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ABSTRACT 

(  Earthquake  engineers  are  often  faced  with  the  problem  of 
determining  the  temporal  and  spatial  variation  of  near-surface  seismic 
motions  in  a  site.  This  type  of  information  is  needed  for  the 
evaluation  of  soil-structure  interaction  effects,  liquefaction 
potential  and  the  effects  of  local  site  conditions  on  surface  motions. 

Actual  ground  motions  are  due  to  a  complicated  system  of  body 
waves  and  surface  waves.  However,  it  is  usually  assumed  that 
near-surface  motions  consist  only  of  vertically  propagating  waves.  In 
order  to  examine  the  validity  of  this  assumption  for  engineering  design 
a  theoretical  investigation  has  been  made  into  the  nature  of 
near-surface  actions  produced  by  horizontally  propagating  waves.  These 
include  inclined  P-,  SV-,  and  SH -waves ,  Rayleigh  waves  and  Love  waves 
in  horizontally  layered  sites  over  a  viscoelastic  half  space. 

The  research  involved  five  phases;  (1)  review  of  current 
knowledge,  (2)  development  of  new  methods  of  site  response  analysis, 

(3)  application  to  site  response  analysis,  (4)  application  to 
soil-structure  interaction  analysis  and  (5)  evaluation  of  the  relative 

v 

importance  of  horizontally  propagating  waves  in  engineering  design. 

The  new  method  of  site  response  analysis  involves  a  finite  element 
type  discretization  of  the  site  in  the  vertical  direction.  According 
to  this  method  the  site  is  subdivided  into  thin  sublayers  and  it  is 
assumed  that  displacements  vary  linearly  between  layer  interfaces.  The 
method  is  essentially  linear  and  works  in  the  frequency  domain. 


Nonlinearities  are  handled  by  an  equivalent  linear  method  according  to 
which  the  stiffness  and  damping  ratio  within  each  layer  are  adjusted 
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iteratively  to  be  compatible  with  the  strains  developed  in  the  layer. 
Transient  motions  are  handled  by  Fourier  techniques. 

The  method  is  essentially  the  same  for  inclined  body  waves  and 
surface  waves.  However,  the  latter  requires  the  solution  of  a  special 
eigenvalue  problem  and  identification  of  the  fundamental  mode.  The 
procedures  have  been  implemented  in  the  two  computer  codes,  SITE  and 
LOVE.  These  codes  can  produce  the  complete  transient  field  of  motion 
from  the  knowledge  of  the  motion  at  one  point  and  the  type  of  wave 
field  producing  the  motion.  Any  specified  combination  of  inclined  body 
waves  and  surface  waves  can  be  considered. 

The  procedure  has  been  applied  to  a  number  of  sites  (rock,  sand, 
and  alluvium)  assuming  different  types  of  wave  fields  and  the  motions 
produced  by  these  fields  are  compared  with  those  produced  by  vertically 
propagating  waves.  The  results  show  that  the  realistic  analysis  of 
incident  body  waves  produce  near-surface  motions  which  vary  with  depth 
in  essentially  the  same  manner  as  those  produced  by  vertically 
propagating  body  waves.  The  motions  produced  by  surface  waves  are 
somewhat  different.  However,  the  study  shows  that  in  soil  sites 
surface  wave  motions  decay  rapidly  in  the  direction  of  wave  propagation 
Within  a  few  hundred  feet  of  the  control  motion  all  components  of 
frequencies  higher  than  1  Hz  are  reduced  to  insignificant  amplitudes. 
The  same  phenomenon  occurs  in  rock  sites  but  at  a  much  slower  rate.  It 
is  therefore  questionable  whether  high  frequency  surface  waves  are 
important  for  engineering  design. 

The  study  of  soil-structure  interaction  effects  show  that  only  the 


free-field  motions  within  the  body  of  soil  replaced  by  the  structure 
are  of  importance  for  design.  Thus,  with  a  specified  surface  control 


notion,  the  spatial  variation  of  free-field  ground  motions  need  to  be 
determined  only  within  a  shallow  zone  near  the  surface. 

Examples  of  soil -structure  interaction  analyses  are  provided  for  a 
structure  on  rock,  a  structure  on  sand,  and  a  large  retaining  wall  on 
an  alluvial  site.  The  results  show  that  for  realistic  wave  fields  the 
motions  of  structures  on  soil  site  depends  only  slightly  on  the  type  of 
wave  field  assumed.  On  rock  sites  surface  waves  may  produce  somewhat 
larger  motions  than  vertically  propagating  body  waves. 

The  study  also  shows  that,  while  the  motions  produced  in 
structures  by  different  types  of  wave  fields  may  not  be  too  different, 
the  dynamic  earth  pressures  on  embedded  structures  depend  strongly  on 
the  nature  of  the  seismic  environment.  In  particular  Rayleigh  waves 
may  produce  larger  dynamic  earth  pressures  than  vertically  propagating 
shear  waves. 

The  final  conclusions  of  the  study  are  that  the  current  assumption 
of  vertically  propagating  waves  is  probably  sufficient  for  many 
practical  purposes.  However,  surface  waves  may  be  important  in  rock 
sites  and  in  the  determination  of  forces  acting  on  structures. 
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CHAPTER  1 
INTRODUCTION 

Earthquake  engineers  are  often  faced  with  the  problen  of 
determining  the  spatial  and  temporal  variation  of  seismic  motions  in  a 
soil  profile  from  a  motion  specified  at  a  single  point.  Solutions  to 
such  problems,  which  are  known  as  site  response  problems,  are  necessary 
for  liquefaction  and  soil-structure  interaction  analyses. 

1.1  Current  Methods  of  Site  Response  Analysis 

Current  engineering  analyses  of  site  response  usually  involve 
three  basic  assumptions: 

•  Ground  motions  developed  near  the  surface  of  a  soil  deposit  may 
be  attribute  only  to  the  vertical  propagation  of  shear  waves, 
Kanai  (1950,  1952). 

e  The  ground  surface,  the  interfaces  between  layers,  and  the 
bedrock  are  essentially  horizontal. 

e  The  material  in  each  layer  is  homogeneous  and  linearly  elastic 
or  viscoelastic. 

Using  these  assumptions,  many  researchers  have  developed 
computational  site  models  and  methods  of  analysis  for  site  response 
problems,  including  Xdriss  and  Seed  (1967),  Tsai  (1969),  Roesset  and 
Whitman  (1969),  and  Schnabel,  Lysmer,  and  Seed  (1972).  The  first  two 
assumptions  above  were  found  to  be  quite  reasonable  for  many  sites 
involving  sedimentary  deposits  with  horizontal  layering.  The  third 
assumption  of  linearity  might  be  inappropriate  for  strong  seismic 
motion.  However,  the  nonlinear  behavior  of  soil  can  be  practically 
approximated  by  the  equivalent  linear  method  proposed  by  Seed  and 


Idriss  (1969). 
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In  general,  the  computational  methods  for  site  response  problems 
can  be  separated  into  continuun  or  discrete  methods.  Continues  methods 
involve  either  the  analytical  solution  of  boundary  value  problems 
directly  frcm  differential  equations,  the  method  of  characteristics,  or 
the  finite-difference  method.  These  solutions  can  be  obtained  either 
in  the  frequency  or  the  time  danain.  Discrete  methods  use  lumped-mass 
or  consistent-mass  finite  element  formulations ,  which  give  good  results 
if  each  layer  in  the  model  is  thin  enough  to  transmit  the  shortest 
wavelength  involved.  Current  discrete  methods  include  the  oomplex 
response  method,  modal  analysis,  direct  integration,  and  the  method  of 
characteristics.  Each  method  is  briefly  reviewed  below. 

Complex  response  analysis  (linear  frequency  domain  analysis)  can 
conveniently  account  for  material  damping  through  the  introduction  of 
complex  moduli  into  the  equations  of  motion.  This  method  can 
incorporate  equivalent  linear  techniques  to  approximate  nonlinear  soil 
behavior.  Furthermore,  Fast  Fourier  Transform  and  interpolation 
techniques  in  the  frequency  domain  make  this  method  effective  and 
economical. 

Modal  analysis  can  be  performed  on  a  lumped-mass  model  of  a  shear 
beam  representing  the  soil  profile.  Modal  frequencies  and  mode  shapes 
may  be  obtained  from  the  geometry  and  mass  distribution  of  the  system. 
The  response  in  each  mode  may  then  be  determined,  and  the  total 
response  is  obtained  by  superposition.  This  technique,  however,  can 
not  properly  account  for  the  spatial  variation  of  damping  within  the 
soil  mass  or  for  radiation  damping. 

Direct  integration  (step-by-step  time  domain  analysis)  may  also  be 
used.  Experience  has  shown  that  for  a  large  time  step  this  method 
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encounters  stability  and  damping  problems  related  to  the  numerical 
algorithms.  Although  mealier  time  steps  will  overcame  same  of  these 
problems,  computational  costs  increase  dramatically.  Problems  may  also 
be  encountered  in  constructing  the  proper  damping  matrix  from  given 
material  properties.  Besides,  unreliable  responses  might  be  obtained 
in  the  high-frequency  ranges. 

The  method  of  characteristics  is  a  mathematical  technique  for 
transforming  partial  differential  equations  into  ordinary  differential 
equations  that  are  then  solved  by  seme  suitable  technique.  The  method 
is  effective  for  linear  analysis  but  encounters  computational  problems 
far  nonlinear  analysis. 

Recently,  efforts  have  been  directed  toward  the  development  of 
better  nonlinear  analysis  methods.  Nonlinear  total  stress  analyses  of 
site  response  problems  using  vertically  propagating  shear  waves  were 
studied  by  Streeter  et  al.  (1974),  Constantopoulos  (1973),  Papadakis 
(1973),  Joyner  and  Chen  (1975),  Martin  (1975),  Idriss  et  al.  (1976), 
and  Taylor  and  Larkin  (1978) .  Several  methods  of  nonlinear  effective 
stress  analysis  have  been  proposed  by  Finn  et  al.  (1977),  Liou  et  al. 
(1977) ,  and  Ghaboussi  and  Dikmen  (1978) .  These  methods  are  important 
for  the  study  of  liquefaction  during  earthquakes. 

1.2  The  Seismic  Environment 

The  above  methods  assume  a  simple  seismic  environment  consisting 
of  vertically  propagating  waves.  However,  as  shown  in  Fig.  1.1  the 
motions  generated  by  a  source  in  the  Earth's  crust  are  composed  of  many 
wave  types.  The  basic  wave  types  are  called  shear  waves  (S-waves)  and 
oompressional  waves  (P-waves) .  When  the  (perhaps  curved)  wave  fronts 
of  these  waves  impinge  on  the  ground  surface  or  layer  interfaces 


surface  waves  may  be  generated.  These  include  Rayleigh  waves  (R-waves) 
and  Love  waves  (L-waves) . 

The  different  wave  types  can  be  classified  as  shown  in  Fig.  1.2. 
S-waves  involve  motions  perpendicular  to  the  direction  of  wave 
propagation.  S-wave  motions  in  the  vertical  plane  are  called  SV-waves. 
Horizontal  S-wave  motions  are  called  SH-waves.  P-waves  involve  motions 
in  the  direction  of  wave  propagation.  Rayleigh  waves  involve 
horizontally  propagating  elliptical  motions  in  the  vertical  plane  and 
Love  waves  consist  of  horizontal  motions  perpendicular  to  the 
horizontal  direction  of  wave  propagating. 

All  of  the  above  wave  types  can  propagate  independently.  However, 
at  layer  interfaces  and  other  inhcmogeneties  refraction  or  reflection 
may  occur  which  not  only  may  change  the  direction  of  wave  propagation 
but  which  may  convert  one  wave  type  into  another  (mode  conversion).  As 
a  result  actual  seismic  environments  are  much  more  oomplicated  than  the 
vertically  propagating  wave  field  assumed  in  current  engineering 
analyses. 

The  main  purpose  of  the  research  described  herein  was  to 
investigate  the  possibility  of  developing  methods  of  site  response 
analysis  which  can  consider  more  realistic,  and  thus  more  oomplicated, 
seismic  environments  than  that  described  in  Section  1.1.  Specifically, 
the  assumption  of  vertical  wave  propagation  will  be  dropped. 

1.3  Horizontally  Propagating  Waves 

Five  types  of  horizontally  propagating  wave  fields  in  horizontally 
layered  soil  and  rock  systems  will  be  investigated: 

•  Inclined  P-waves 


•  Inclined  SV-waves 


SHEAR  WAVES  AND  SHEAR 


Fig.  1.1  Idealized  Relation  Among  Earthquake  Source,  Wave  Paths 
and  Site  (after  Taal,  1969:  Nalr,  1975) 
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Fig.  1.2  Seismic  Waves  -  Types,  Characteristics  and  Relationships 


•  Inclined  SH- waves 


•  Rayleigh  waves 

•  Love  waves 

Also,  combinations  of  such  wave  systems  will  be  considered. 

The  relative  contributions  of  the  different  wave  types  to  the 
total  ground  motion  and  the  corresponding  arrival  times  for  each  type 
of  wave  depend  on  the  epicentral  distance  to  the  site,  the  focal  depth 
of  the  source,  and  the  phenomena  of  multiple  reflection,  refraction,  and 
dispersion  along  the  various  paths.  Realistically,  it  must  be  assumed 
that  all  observed  seismograms  contain  some  oomponents  of  all  of  the 
above  motions.  However,  the  exact  composition  is  not  and  probably 
never  will  be  known  since  the  exact  properties  of  the  source  and  the 
physical  and  geometric  details  of  the  geology  cannot  be  determined. 

There  is  today  considerable  observational  evidence  that  all  of  the 
above  wave  types  exist.  The  existence  of  inclined  body  waves  has  been 
confirmed  by  numerous  investigators,  e.g.  Suzuki  (1932)  who  reported  a 
mean  angle  of  incidence  of  about  4°  for  the  initial  motion  in  about 
fifty  records  from  Hongo  and  Mitika,  Japan. 

The  evidence  for  the  existence  of  surface  waves  is  even  stranger. 
More  recent  observations  include;  Shima  (1970),  Trifunac  (1971.,  Bolt 
(1972),  Anderson  (1974),  Hanks  (1975),  Toki  (1977),  Swanger  and  Boore 
(1978).  However,  the  evidence  seems  to  be  limited  to  frequencies  below 
2  Hz.  For  example,  the  surface  waves  observed  by  Swanger  and  Boore 
(1978)  in  the  records  from  the  1968  El  Centro  earthquake  occurred  in 
the  range  0.1  -  1.0  Hz  and  the  surface  waves  detected  in  the  1971  San 
Fernando  earthquake  were  in  the  range  0.1  -  2.0  Hz,  Toki  (1977).  The 
existence  of  significant  surface  wave  oomponents  above  2  Hz,  which  is 


j 
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the  range  of  most  interest  to  engineers,  has  therefore  not  been 
confirmed. 


In  addition  to  the  above  evidence  there  is  an  overwhelming  amount 
of  data  which  indicate  phase  differences  between  motions  observed  at 
closely  spaced  points,  e.g.  Yamahara  (1970)  who  observed  significant 
phase  differences  between  stations  spaced  only  about  100  feet  apart. 

In  most  of  these  cases  the  wave  type  was  not  identified.  Nevertheless, 
the  evidence  confirms  the  existence  of  horizontally  propagating  waves. 

Although  a  few  suggestions  have  been  made,  Nair  and  Emery  (1975), 
Liang  and  Duke  (1978),  as  to  the  relative  content  of  different  wave 
types,  the  literature  has  a  dearth  of  data  on  this  topic  in  the 
frequency  range  of  interest  to  engineers. 

Several  researchers  have  developed  theories  for  the  response  of  a 
horizontally  layered  site  to  plane  harmonic  body  waves  arriving  at  a 
specified  incident  angle  from  an  underlying  half  space.  Thomson  (1950) 
and  Haskell  (1960,  1962)  developed  a  matrix  formulation  for  computing 
transmission  coefficients  in  a  layered  continuum.  Hannon  (1964) 
extended  Haskell's  formula  ion  to  study  transient  incident  P-waves. 
Silva  (1976)  extended  the  Thomson-Haskell  method  to  include  damping  in 
soil  layers. 

The  response  of  a  horizontally  layered  site  to  harmonic  surface 
waves  has  been  studied  by  Sezawa  and  Kanai  (1935),  Haskell  (1953),  and 
Ewing  et  al.  (1957),  Mooney  and  Bolt  (1966).  Recently,  Bocheva  (1977) 
extended  the  method  to  study  surface  wave  amplification  factors. 

Swanger  and  Boore  (1978)  simulated  strong  motion  displacement  using 
surface  wave  modal  superposition.  Lyamer  (1969a),  Waas  (1972),  Lyaner 
and  Naas  (1972),  and  Lysmer  and  Drake  (1972)  applied  the  finite  element 
method  to  problems  involving  Rayleigh  and  Love  waves. 
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Most  of  the  above  studies  are  restricted  to  linear  analysis  of 
harmonic  motion  of  a  single  type  of  wave  field. 

1.4  Purpose  and  Scope 

The  main  purpose  of  the  research  described  herein  was: 

•  To  develop  analytical  methods  for  site  response  analysis  of 
horizontally  layered  sites  excited  by  horizontally  propagating 
seismic  motions  consisting  of  surface  waves  and  inclined  body 
waves . 

The  main  emphasis  is  on  engineering  applications.  This  means  that 
the  geometric  dimensions  of  the  model  considered  are  smaller  than  (and 
the  frequency  range  higher  than)  those  usually  considered  by 
seismologists  (100  ft  vs.  1  km,  1-20  Hz  vs  0.1  -  1  Hz).  Also,  while 
the  seismologist's  problem  usually  is  to  determine  motions  from 
estimated  source  parameters  or  source  parameters  from  observed  motions, 
the  engineering  site  response  problem  involves  determining  the  spatial 
and  temporal  variations  of  transient  motions  within  a  limited  distance 
from  a  specified  motion  (control  motion)  at  or  near  the  ground 
surface.  This  process  is  called  deconvolution  in  the  engineering 
profession. 

The  research  involved  the  following  items: 

•  Review  of  existing  methods  and  available  data  on  dynamic 
material  properties. 

•  Development  of  a  finite  element  method  for  transient  site 
response  analysis  of  problems  involving  inclined  body  waves  in 
a  profile  consisting  of  soil  layers  over  a  viscoelastic  half 


ai 
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space. 
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•  Development  of  a  finite  element  method  for  transient  site 
response  analysis  of  problems  involving  Rayleigh  waves  and  Love 
waves  in  a  profile  consisting  of  soil  layers  over  a 
viscoelastic  half  space. 

•  Development  of  practical  computer  codes  (SITE  and  LOVE)  to 
implement  the  above  methods. 

•  Application  of  the  above  methods  to  realistic  site  response 
problems. 

•  Application  of  site  response  solutions  to  soil-structure 
interaction  problems. 

The  presentation  of  the  research  is  organized  as  follows: 

The  dynamic  properties  of  soils  and  rock  are  described  in 
Chapter  2.  The  emphasis  is  on  a  material  description  which  is  suited 
for  the  analytical  procedures  employed  in  later  chapters. 

Inclined  body  waves  are  discussed  in  Chapter  3  and  surface  waves 
in  Chapter  4.  The  treatment  in  these  chapters  involves  only  harmonic 
waves.  However,  these  chapters  contain  most  of  the  theoretical 
developments  for  the  finite  element  codes  SITE  and  LOVE. 

The  transition  to  transient  cases  through  Fourier  techniques  is 
made  in  Chapter  5  which  also  contains  several  case  studies  of  site 
response  analysis. 

In  Chapter  6  the  application  of  site  response  solutions  to 
soil-structure  interaction  problems  is  discussed  and  a  number  of  case 
studies  are  presented. 

The  results  of  the  research  are  summarized  in  Chapter  7.  As  will 
be  shown  in  that  chapter  the  research  lead  not  only  to  the  development 
of  a  unified  theory  for  inclined  body  waves  and  surface  waves  in 


layered  systems  and  two  associated  ooaputer  programs  but  to  a  number  of 
significant  conclusions  regarding  the  importance  of  considering 
horizontally  propagating  waves  in  design  and,  perhaps  surprisingly,  the 
likely  contribution  of  surface  waves  to  near  surface  seismic  motions. 


11 


CHAPTER  2 

DYNAMIC  MATERIAL  PROPERTIES 


2.1  Introduction 

The  stress-strain  characteristics  of  soils  are  strongly  nonlinear 
and  may  significantly  influence  the  dynamic  response  of  a  site  subjected 
to  a  strong  earthquake.  A  good  site  response  analysis  must  therefore 
consider  these  nonlinear  effects. 

Details  of  the  dynamic  stress-strain  behavior  of  soils  have  recently 
been  reported  in  state-of-the-art  papers  by  Hardin  (1978)  and  Yoshimi 
et  al.  (1977) .  It  is  clear  from  these  reports  that  the  transient 
stress-strain  behavior  of  soils  is  extremely  complicated  and  that  this 
behavior  can  not  as  yet  be  fully  described  by  constitutive  laws.  Most 
of  the  data  and  models  available  refer  to  cyclic  behavior  of  soils 
subjected  to  constant  strain  amplitudes.  Typical  stress-strain 
relationships  of  soils  subjected  to  symmetric  cyclic  loading  conditions 
are  curvilinear  as  shown  in  Pig.  2.1. 

In  choosing  dynamic  soil  properties  for  site  response  analysis,  one 
should  realize  that  such  problems  can  only  be  solved  by  making  certain 
assumptions  about  the  nature  of  the  wave  fields  involved.  Except  for 
the  special  case  of  vertically  propagating  waves,  which  is  not  the 
major  topic  of  this  dissertation,  appropriate  wave  fields  can  only  be 
constructed  for  linear  layered  systems.  Therefore,  it  is  essential  to 
choose  representative  linear  dynamic  properties  for  the  actual  analysis. 
As  will  be  shown,  such  properties  can  be  determined  from  the  available 
data,  and  an  approximation  to  nonlinear  analysis  can  be  achieved  by  the 
equivalent  linear  method,  which  is  discussed  at  the  end  of  this  chapter. 


... . 
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2.2  Measurement  of  Dynamic  Properties 

Considerable  effort  has  baan  directed  towards  tha  datarmination  of 
soil  propartias  in  racant  years.  A  complete  review  of  the  measurement 
of  dynamic  properties  was  given  in  a  state-of-the-art  paper  by  woods 
(1978).  The  most  commonly  used  test  procedures  are  described  below. 
Determination  of  Hysteresis  Loops 

Hysteretic  stress-strain  relationships  of  the  type  shown  in 
Fig.  2.1  can  be  determined  by  cyclic  triaxial  compression  tests,  cyclic 
simple  shear  tests,  or  cyclic  torsional  shear  tests.  These  tests  are 
applicable  in  the  amplitude  ranges  shown  in  Fig.  2.2a  and  are  usually 
performed  in  the  frequency  range  of  1-20  Hz.  Test  results  have 
indicated  that  the  shape  of  the  hysteresis  loops  is  virtually 
independent  of  frequency.  From  these  loops  the  effective  dynamic 
moduli  and  fractions  of  critical  damping  can  be  determined.  The 
modulus  is  the  slope  of  the  secant  betweeen  the  ends  of  the  loop,  and 
the  damping  is  proportional  to  the  aspect  ratio  of  the  loop,  i.e.,  the 
ratio  between  the  average  width  and  the  length  of  the  loop. 

Resonance  Column  Tests 

The  dynamic  moduli  can  also  be  determined  from  longitudinal  or 
tortional  resonance  tests,  in  which  a  column  of  soil  is  excited  at 
different  frequencies  to  determine  the  natural  frequencies  from  which 
the  moduli  can  be  computed.  The  damping  ratio  can  be  estimated  from 
the  height  of  the  resonance  peaks  or  by  measurement  of  the  phase 
difference  between  the  displacement  of  the  specimen  and  the  exciting 
force.  These  tests  are  usually  performed  at  frequencies  in  the  range 
of  20-260  Hz  and  are  applicable  in  the  strain  amplitude  ranges  shown  in 
Fig.  2.2a. 


Hysteretic  Stress-Strain  Relationships 

at  Different  Strain  Amplitudes  Pig-  2.2  Shearing  Strain  Amplitude  Capabilities  for 

Laboratory  tests  and  Field  Techniques 
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Fr»«  Vibration  Tests 

In  these  tests  cylindrical  soil  samples  are  set  into  longitudinal 
or  torsional  vibrations,  the  power  is  switched  off,  and  the  decay  of 
the  amplitudes  is  measured  to  determine  the  logarithmic  decrement,  from 
which  the  damping  ratio  can  be  computed.  These  tests  can  be  conducted 
using  a  resonant  column  apparatus  with  either  solid  or  hollow  samples, 
and  good  results  can  be  obtained  only  at  relatively  low  to  moderately 
high  strain  levels. 

Field  Measurement  of  Wave  Velocities 

In-situ  tests  are  conducted  to  determine  the  velocities  of 
propagation  of  P- ,  S-  and  Rayleigh  waves.  The  most  usual  types  of 
field  tests  are 

a.  Geophysical  tests:  seismic  refraction,  seismic  cross-hole,  and 
seismic  down-hole  methods. 

b.  Surface  vibration  tests:  surface  wave  and  resonant  footing 
techniques. 

c.  Other  field  techniques:  frequency  domain  measurements; 
cylindrical  in-situ  tests,  see  Woods  (1978) . 

In  general,  these  tests  give  soil  moduli  for  low  strain  levels. 
In-situ  tests  are  inadequate  to  determine  the  volumetric  characteristics 
of  saturated  soil  because  the  measured  P-wave  velocities  are  greatly 
affected  by  the  presence  of  water.  The  customary  procedure  in  this 
case  is  to  conduct  laboratory  measurements  of  Poisson's  ratio  or  bulk 
modulus.  Such  laboratory  tests  are  usually  performed  statically  because 
the  measurement  of  lateral  deformations  and  volumetric  strains  under 
dynamic  conditions  is  not  practical.  Shannon  and  Wilson  (1971). 

The  ranges  of  shear  strain  amplitudes  over  which  the  field 


techniques  are  applicable  are  shown  in  Fig.  2.2b. 
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2.3  Suaaarv  of  Available  Data 

As  can  be  seen  from  Fig.  2.1,  the  effective  dynamic  moduli  of  soils 
subjected  to  cyclic  excitation  will  usually  decrease  with  the  strain 
amplitude,  and  the  damping  ratio  will  increase  as  the  strain  amplitude  is 
increased.  Therefore,  it  is  customary  to  present  the  data  in  the  form  of 
curves  which  show  the  variation  of  modulus  and  damping  with  strain 
amplitude,  see  Figs.  2. 3-2. 9. 

Most  available  data  only  consider  the  variation  of  shear  modulus  and 

damping  ratio  with  shear  strain  amplitude.  In  principle,  two  moduli  and 

two  damping  ratios  should  be  considered — one  set  corresponding  to  S-waves 

(shear  modulus,  G,  and  the  corresponding  damping  ratio,  Bg) ,  and  a 

second  set  corresponding  to  P- waves  (constrained  modulus,  M,  and  damping 

ratio,  3  ).  Also,  the  variation  of  these  parameters  with  the  amplitude 
P 

of  normal  strain  should,  in  principle,  be  considered.  However,  normal 
strains  are  usually  considered  to  have  only  small  effects  on  the  dynamic 
properties  and  are  neglected.  Similarly,  the  two  damping  ratios  Bg  and 
Bp  are  often  assumed  to  be  the  same,  although  data  by  McDonal  (1958) 
and  Eisenburg  (1972)  indicate  that  the  damping  ratio  for  S-waves  is 
considerably  higher  than  that  for  P-waves  (see  below).  A  detailed  study 
of  the  factors  influencing  the  shear  moduli  and  damping  values  of  soils 
has  been  carried  out  by  Hardin  and  Drnevich  (1972).  Seed  and  Idriss 
(1970)  have  proposed  simplified  practical  relations  which  will  be  used  in 
this  study.  Sane  relationships  for  rock  material  were  proposed  by 
Schnabel  (1973). 

Cohesionless  Soils 

Seed  and  Idriss  (1970)  have  shown  that  the  dynamic  shear  modulus  of 
cohesionless  soil  can  be  expressed  by: 

G  ■  1000  K,  (O  )“  F 


(2.1) 
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where  is  a  parameter  that  depends  on  the  relative  density  and  the 
shear  strain  amplitude.  Also,  0  is  the  mean  effective  stress,  which 

!D 

equals  (0v  +  2XQ  oh)/3,  where  a v  and  0h  are  the  vertical  and 
horizontal  effective  pressures,  respectively,  and  XQ  is  the  at-rest 
earth-pressure  coefficient.  The  exponent,  a,  has  been  found  to  vary  from 
0.3  to  0.8  (Idriss  and  Seed,  1968;  Carriveau,  1970;  Drnevich  et  ad., 

1966;  and  Silver  and  Seed,  1969).  Seed  and  Idriss  (1970)  proposed  the 
use  of  a  ■  0.5  in  the  above  expression.  The  term  F  is  a  coefficient 
accounting  for  grain  size  and  shape  variation.  It  ranges  from  0.6  for 
silt  to  2  for  gravel. 

The  estimated  average  value  for  the  combined  effects  of  K2  and  F 

-4 

is  about  61  at  low  shear  strain  levels  (10  percent)  for  a  wide  range 

of  sandy  soils  at  75%  relative  density.  The  attenuations  of  the  shear 
modulus  with  increasing  strain  for  sands  of  different  densities  are  shown 
in  Fig.  2.3.  The  curves  shown  in  this  figure  may  be  normalized  to  a 
single  attenuation  curve  as  shown  in  Fig.  2.4. 

Measured  values  of  the  damping  ratio  for  cohesionless  soils  and 
average  values  proposed  by  Seed  and  Idriss  (1970)  are  shown  in  Fig.  2.5. 
The  average  curve  shown  is  adequate  for  most  cohesionless  soils  up  to  a 
confining  pressure  of  2500  psf.  The  damping  ratio  will  be  affected  by 
overburden  pressure,  relative  density,  degree  of  saturation  and  the 
number  of  loading  cycles.  The  effects  of  the  number  of  loading  cycles 
and  relative  density  are  minor.  It  has  been  found  that  the  damping 
decreases  with  increasing  effective  overburden  pressure  and  increases 
with  increasing  degree  of  saturation.  Schnabel  (1973)  suggested  that  the 
variation  of  the  damping  factor  with  effective  overburden  pressure,  a  , 
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Variation  of  Shear  Modulus  with  Shear  Strain  for  Sands 
(after  Seed  and  Idriss,  197) 
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may  be  expressed  by 
1/4 


6 


2500 
a 


\  6.y 
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where  Fs  is  0.6  for  dry  sand  or  1.0  for  saturated  sand. 
Cohesive  Soil 


The  shear  modulus  of  saturated  clays  has  been  found  to  be 


(2.2) 


essentially  proportional  to  the  undrained  shear  strength  Su»  as  follows: 


G*KSu  (2.3) 

where  K  is  a  coefficient  depending  mainly  on  the  shear  strain  amplitude. 
Average  values  of  the  coefficient  K  obtained  from  various  sources  are 
shown  in  Fig.  2.6,  Seed  and  Idriss  (1970).  Some  recent  results  by  Stokoe 
and  Lodde  (1978)  are  shown  in  Fig.  2.7.  The  shear  modulus  of  clay  can  be 
determined  in  the  field  for  low  strain  levels  and  in  the  laboratory  for 
high  strain  ranges.  Sample  disturbance  will  significantly  affect  the 
shear  moduli  obtained  from  laboratory  tests.  Hence,  laboratory  data  must 
be  corrected  for  sample  disturbance.  Correction  factors  can  be  obtained 
by  comparison  of  low  strain  tests  in  both  field  and  laboratory. 

The  damping  ratio  of  clays  is  affected  by  shear  strain  amplitude, 
effective  overburden  pressure,  void  ratio,  number  of  loading  cycles  and 
water  content.  Most  available  data  cover  only  the  effect  of  shear  strain 
amplitude.  Measured  values  of  the  damping  ratio  and  proposed  average 
values  for  saturated  clay  obtained  from  different  sources  of  data  are 
shown  in  Fig.  2.8. 


Rock 

Values  of  shear  modulus  for  rock  are  most  often  obtained  from 
seismic  investigations,  which  yield  values  only  at  low  strain  levels. 
Very  few  data  are  available  for  strain  dependence,  but  it  seems  likely 
that  rock  will  exhibit  some  decrease  of  shear  modulus  with  increasing 
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Fig.  2.7  Norulltad  Shear  Modulus  with  Shear  Strain  for 
Saturated  Clays  (froa  Stokoe  and  Lodde,  1978) 


Fig.  2.8  Deaping  Ratios  for  Clays  (Seed  and  Idriss) 
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shear  strain.  The  curves  shown  in  Fig.  2.9  were  proposed  by  Schnabel 
(1973)  for  site  response  analysis  of  sedimentary  rock  layers.  Schnabel 
also  presented  considerable  data  for  1cm  strain  properties  of  different 
rock  types. 

A  literature  survey  by  Knopoff  (1964)  provides  values  of  damping 
ratios  for  various  rock  types.  These  data  were  obtained  from  laboratory 
tests  and  indicate  that  the  damping  ratio  for  rock  varies  from  0  to  1.4 
percent.  A  literature  survey  by  Jackson  and  Anderson  (1970)  of  data 
obtained  from  in-situ  measurement  in  the  shallow  crust  of  the  earth 

-2 

indicates  values  for  the  S-wave  damping  ratio  ranging  from  2.5  x  10  % 

to  0.5%,  and  values  for  the  P-wave  damping  ratio  ranging  from  0.01%  to 
1.5%.  For  surface  earth  materials,  Knopoff 's  survey  (1964)  shows  values 
of  the  P-wave  damping  ratio  ranging  from  1%  for  magnetite  hemotite  to 
7.2%  for  Pottsville  sandstone  and  values  of  the  S-wave  damping  ratio  of 
about  5%  for  Pierre  Shale. 

Damping  for  shear  waves  is  higher  than  for  P-waves  by  a  factor  of 

1.8  to  2.6  (McDonal,  1958;  Eisenberg,  1972).  No  data  are  available 

regarding  the  variation  of  damping  with  strain  in  rock,  although  some 

increase  in  damping  ratios  with  increasing  strain  is  to  be  expected.  The 

strain  dependent  damping  curve  shown  in  Fig.  2.9  was  proposed  by  Schnabel 

(1973)  for  sedimentary  rock  layers  with  shear  wave  velocities  in  the 

« 

range  2000-4000  fps  at  100-3000  ft  depth.  Be  also  proposed  the  following 
method  for  adjusting  the  damping  values  of  rock  having  other  shear  wave 
velocities. 

a.  All  materials  with  V^<3000  fps  are  treated  as  soil  with  the 
same  nonlinear  property  characteristics  as  described  for  soils. 

b.  All  materials  with  Vs>ll,000  fps  may  be  treated  as  linear 
elastic  materials. 
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c.  A  material  with  3000  fps<Vs<n,000  fps  is  treated  as  a 

transition  material  between  soil  and  linear  elastic  rock,  in 
this  case,  linear  interpolation  can  be  employed  between  the  two 
extremes  described  above. 

Dynamic  Poisson's  Ratio 

The  dynamic  Poisson's  ratio  of  soil  during  cyclic  loading  has 
attracted  very  little  attention.  Fig.  2.10  shows  some  laboratory  test 
results  for  the  dynamic  Poisson's  ratio  of  clayey  soil  tested  at 
different  shear  strain  levels,  Hara  (1973) .  It  can  be  seen  that  the 
dynamic  Poisson's  ratio  for  a  soft  clay  and  a  very  stiff  clay  are 
essentially  independent  of  the  shear  strain  levels  and  frequencies  of 
cyclic  loading.  A  statistical  analysis  of  recorded  wave  velocities  of 
various  deposits  obtained  by  seismic  exploration  was  conducted  by  Ohsaki 
and  Iwasaki  (1973).  The  evaluated  dynamic  Poisson's  ratio  versus  shear 
moduli  and  the  total  average  values  are  shown  in  Fig.  2.11.  The  results 
demonstrate  that  the  dynamic  Poisson's  ratio  does  not  change  appreciably 
with  shear  moduli  for  sandy  soils.  Using  the  experimental  results  given 
by  Hara  (1973),  they  concluded  that  the  dynamic  Poisson's  ratio  for 
cohesive  soils  is  almost  constant  (approximately  0.48)  and  that  for 
cohesionless  soil,  the  dynamic  Poisson's  ratio  is  a  function  of  shear 
modulus. 

In  fact,  the  dynamic  Poisson's  ratio  will  significantly  affect 
stress  and  strain  computations  as  well  as  the  wave  propagation 
characteristics  of  P-waves  in  a  soil  deposit  during  dynamic  excitation. 
2.4  Theoretical  Models 


Several  theoretical  constitutive  models  have  been  proposed  which, 
with  cyclic  excitation,  produce  hysteresis  loops  similar  to  those  shown 
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Pig.  2.10  Dynamic  Poisson's  ratio  of  clay 
(after  Hara,  1973) 
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Pig.  2.11  Dynaaic  Poisson's  ratios  ar.i 

shear  soduli  (after  ohsaki,  1973) 


in  Fig.  2.1.  For  the  special  case  of  the  shear  stress-shear  strain 
relationship  the  most  relevant  models  for  soil  dynamics  analysis  are 
those  discussed  below. 

1.  Viscoelastic  Models 

The  general  viscoelastic  solid  produces  an  elliptical  hysteresis 
loop,  the  shape  and  slope  of  which  are  rate  dependent  (i.e.,  the 
effective  modulus  and  damping  ratio  are  frequency  dependent).  However, 
numerous  tests  have  shown  that  the  shapes  of  the  hysteresis  loops  for 
soils  are  essentially  independent  of  frequency  within  the  frequency  range 
of  interest  in  earthquake  engineering.  A  viscoelastic  material  which 
satisfies  this  condition,  as  used  in  this  dissertation,  can  be  defined  by 
the  complex  dynamic  modulus: 

G*  «  G (1  -  28*  +  2iB  Vl  -  b|)  (2.4) 

s  s  s 

%  G (l  +  2iB  )  ;  for  small  8 
s  s 

where  G  is  the  usual  shear  modulus  and  8  is  the  fraction  of  critical 

s 

damping.  With  this  definition  the  stress-strain  law  for  harmonic 
excitation  becomes 
T  »  G*y 

where  T  and  X  are  the  complex  amplitudes  of  stress  and  strain.  A  more 
complex  model  for  two-dimensional  stress  states  will  be  introduced  in 
Chapter  3. 

The  hysteresis  loops  inherent  in  the  above  model  are  independent  of 
the  strain  amplitude.  However,  this  problem  can  be  overcome  by  the 
equivalent  linear  method  discussed  later  in  this  chapter. 

2.  Ramberg-Osqood  Generalized  Model 

A  four-parameter  model  which  can  be  used  for  nonlinear  analysis  was 
proposed  by  Ramberg  and  Osgood  (1943)  and  modified  by  Jennings  (1964) . 
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In  this  model,  the  shear  strain  is  a  function  of  a  given  stress  as 
follows : 


(Y/yy)  -  <T/Tuit)  +a(T/Tult)r  (2,5) 

in  which  r  is  a  positive  constant  greater  than  one,  a  is  a  real  positive 
parameter,  which  is  a  function  of  r,  Yy  is  a  reference  strain,  and 

is  the  ultimate  shear  stress.  The  shear  stress  in  this  equation 
can  not  be  explicitly  represented  by  the  shear  strain  for  a  general  value 
of  r  (r  *  1  gives  a  linear  relationship  between  X  and  T,  and  r  *  •  gives 
an  elasto~plastic  relationship) . 

The  hysteretic  damping,  0,  can  be  evaluated  as  described  by  Jac  en 
(1960) . 


6  ’  8«.x  11  -  t 
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(2.6) 


max 


where  0  ■  (r  -  l)/(irr),  and  G/G  can  be  evaluated  from 

max  max 

G/Gmax  *  1/  [1  +  «(G  r/G„ax  YyJ  (2-7) 

The  model  has  been  used  by  Constantopoulos  and  Christian  (1973)  and 
Streeter  et  al.  (1974)  for  site  response  analysis  and  lately  by  Idriss  et 
al.  (1976)  for  the  gradual  degrading  of  clay  when  subjected  to  cyclic 
loading. 

3.  Hardin-Drnevich  Model 

Hardin  and  Ornevich  (1972b)  proposed  the  following  approximate 
relationship  between  stress  and  strain 


7  *  Gmax  Y/U  ♦  Yh>  (2.8) 

<5  „  can  be  measured  by  resonant  column  or  seismic  techniques  and  y. 
max  n 

is  the  hyperbolic  strain  defined  as: 


Y 


h 


X_ 

Yy 


.  "bY/Y 
jl  ♦  ae 


yi 

I 


(2.9) 


27 


in  which  "a"  and  aba  are  empirical  soil  constants  and  aea  is  the  base  of 
natural  logarithms.  This  model  needs  four  parameters  explicitly.  The 
damping  ratio  is  given  by 


6  ■  6.«  V(1  *  V 


(2.10) 


For  sandy  soil  8  depends  on  the  number  of  cycles  of  loading,  and, 
max 

for  clay  soil,  on  the  loading  frequency  and  stress  state.  A  similar 
model— initially  used  by  Kondner  (1963)  and  lately  also  by  Hardin  and 
Drnevich  (1972b) — does  not  include  the  parameters  a  and  b  and  defines 
Yh  *  y/y  .  This  model  is  called  the  Hyperbolic  Model,  and  only  two 
parameters  are  needed  to  determine  the  stress-strain  relationship. 

4.  Martin-Davidenkov  Model 

Martin  (1975)  modified  the  generalized  Davidenkov  model  by  defining 
a  new  function  for  shear  strain,  and  he  proposed  the  following 
stress-strain  law: 

(2.11) 

where  H(y)  is  given  by 

OT5  _  \  a 

(2.12) 

in  which  A  and  B  are  constant  parameters.  Four  parameters  are  required 
for  this  model.  The  damping  ratio  can  be  evaluated  by 


T  =  l  1  “  H(Y)lY 

max 


H(Y)  -  {(Y/Yy)2B/|l  +  <Y/Yy)2B]]A 


B 


nH(n)  dri}/  {y2 


nH(n)  dn} 


(2.13) 


By  appropriate  choices  of  the  parameters  involved  all  of  the  above 
nonlinear  models  can  be  made  to  fit  approximately  the  strain  dependency 
curves  for  cyclic  loading  published  by  Seed  and  Idriss  (1970),  see 
Fig.  2.12  from  Kagawa  (1978).  The  Martin-Davidenkov  model  provides  the 


closest  fit. 


Parameters  Used  for  Soil  Models 


Hyperbolic  Model  Xy 
Raaberg-Osgood  Model  Yy 
Hardln-Drnevich  Model  Yy 
Martln-Davidenkov  Model  Yy 
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2.5  The  Equivalent  Linear  Method 

The  nonlinear  behavior  of  soil  materials  cannot  be  fully  described 
by  constant  elastic  moduli  and  damping  coefficients.  However,  a  good 
approximation  of  the  effects  of  soil  nonlinearities  on  the  response  can 
be  obtained  by  the  use  of  constant  strain  compatible  moduli  and  damping 
ratios  in  a  sequence  of  linear  analyses.  This  method,  which  is  known  as 
the  equivalent  linear  method  (Seed  and  Idriss,  1969),  can  be  briefly 
described  in  the  following  manner. 

In  a  site  response  analysis  the  equivalent  linear  method  starts  with 
a  linear  analysis  using  estimated  soil  properties  in  each  layer  of  the 
soil  system.  This  analysis  yields  complete  time  histories  of  shear 
strain,  from  which  the  effective  shear  strain  amplitudes  are  calculated 
in  each  layer.  (The  effective  shear  strain  amplitude  is  usually  taken  as 
65%  of  the  maximum  shear  strain  or  as  the  RMS  value  of  the  shear  strain 
time  history).  Using  the  computed  strain  amplitudes,  an  improved  set  of 
soil  moduli  and  damping  ratios  are  obtained  frcm  appropriate  soil  data 
curves  of  the  type  shown  in  Pigs.  2.4-2. 9,  and  a  new  linear  analysis  is 
performed  with  these  properties.  The  process  is  repeated  until  the 
properties  from  two  consecutive  analyses  differ  by  less  than  a  specified 
tolerance,  say  5  percent.  This  will  usually  require  fewer  than  5 
iterations.  The  results  of  the  last  iteration  are  taken  as  the  final 
solution  to  approximate  a  true  nonlinear  solution.  This  technique  has 
been  widely  used  in  practice  because  it  is  an  efficient  method  and  is 
easy  to  implement  in  a  computer  program. 

The  linear  equivalent  method  can  also  be  used  for  two-dimensional 
analysis  by  the  finite  element  method,  Idriss  et  al.  (1973)  and  Lysmer 
et  al.  (1975).  In  such  analyses,  strain  compatible  properties  are 
determined  by  iteration  for  each  soil  element. 
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The  equivalent  linear  technique  has  been  shown  to  give  surprisingly 
good  approximations  to  motions  computed  by  truly  nonlinear  techniques, 
Martin  (1975),  although  cases  have  been  found  where  significant 
differences  occured,  Martin  (1975)  and  Finn  et  al.  (1978).  Good 
agreement  has  also  been  found  between  observed  ground  motions  and  motions 
computed  by  equivalent  linear  methods,  Idriss  and  Seed  (1968),  Schnabel 
and  Seed  (1971),  and  Valera  et  al.  (1977). 

2.6  Combined  Loading  Effect  on  Strain-Dependent  Properties 

The  above  discussion  of  shear  modulus  and  S-wave  damping  is 
basically  for  the  case  of  a  simplified  one-dimensional  S-wave  analysis. 
These  strain  dependent  properties  can  be  easily  obtained  from  available 
laboratory  tests  such  as  the  resonant  column  test  and  the  strain  control 
triaxial  cyclic  test.  However,  during  an  earthquake,  the  soil  is  excited 
simultaneously  by  all  kinds  of  seismic  waves  travelling  in  all 
directions.  An  element  of  soil  will  be  subjected  to  combined  shear  and 
compressive  strains.  How  these  combined  excitations  affect  the  modulus 
and  damping  characteristics  is  still  not  clear. 

A  research  program  to  study  these  combined  effects  in  different 
soils  is  currently  in  progress  at  the  University  of  California  at 
Berkeley.  The  research  conprises  different  studies  in  order  that  a  wide 
range  of  loading  conditions  and  strain  amplitudes  might  be  explored.  For 
the  harmonic  simultaneous  loading  condition,  two  components  of  loading 
can  be  either  out  of  phase  (Rayleigh  wave  type  excitation)  or  in  phase 
excitation  (body  waves  at  small  angles  of  incidence).  At  low  strains, 
studies  involve  the  cyclic  excitation  of  a  triaxial  soil  specimen  in  a 
special  resonant  column  device  capable  of  simultaneous  compression  and 
torsion  excitation.  At  high  strain  ranges,  studies  involve  the  cyclic 
loading  of  a  hollow  cylindrical  specimen  of  soil  with  a  special  testing 
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machine  capable  of  reproducing  earthquake  level  strains.  While  the  study 
is  still  continuing,  the  effect  on  both  the  shear  and  constrained  moduli 
and  damping  appears  to  be  significant  in  seme  high  strain  ranges  and  same 
particular  loading  conditions,  Griffin  (1979). 

Because  of  a  lack  of  experimental  results  on  strain  dependent 
properties  for  simultaneous  loading  conditions,  the  site  response 
analysis  for  R-wave  excitation  is  still  restricted  to  linear  analysis. 
However,  the  developed  computer  program  SITE  can  handle  the  approximate 
nonlinear  analysis  for  R-wave  excitation  whenever  the  strain  dependent 
properties  are  available. 

In  site  response  analyses  with  R-waves  or  inclined  body  waves,  the 
strain- dependent  property  curves  for  both  the  shear  and  the  constrained 
moduli  as  well  as  both  the  S-and  P-wave  damping  ratios  will  be  used 

simultaneously  for  the  iteration  process.  The  complex  constrained 

* 

modulus,  M  ,  is  defined  in  a  similar  way  as  the  complex  shear  modulus 
defined  in  section  2.4s 

M*  -  M(1  -  28  2  +  2iB  Vl  -  B  2)  (2.14) 

P  P  P 

where  M  is  the  real  constrained  modulus  and  Bp  is  the  damping  ratio 
due  to  P-waves.  However,  at  this  stage,  because  of  a  lack  of  data  on 
strain-dependent  P-wave  properties,  one  can  assume  a  constant  real 
Poisson's  ratio  and  follow  the  conventional  iteration  procedures  by 
iterating  on  shear  modulus  and  S-wave  damping.  This  approach  implies 
that  the  analysis  is  using  the  same  rate  of  attenuation  on  both  the  shear 
modulus  and  constrained  modulus  and  also  using  the  same  value  of  damping 
for  S-  and  P-waves.  On  the  other  hand,  as  an  extreme  case  one  can  assume 
a  constant  constrained  modulus  and  a  constant  P-wave  damping  together 
with  strain-dependent  S-wave  properties.  Analyses  iterating  on  these 


sets  of  property  curves  will  result  in  a  complex  Poisson's  ratio.  The 
Poisson's  ratio  will  be  strain  dependent  and  tend  to  be  larger  than  that 
obtained  when  both  6*  and  M*  are  assuaed  to  be  shear  strain  dependent. 


CHAPTER  3 


INCLINED  BODY  HAVES 


3.1  Introduction 

In  this  chapter  a  method  is  developed  for  the  evaluation  of  the 
seismic  response  of  a  horizontally  layered  site  due  to  a  system  of  plane 
incident  body  waves.  These  waves  arrive  at  an  oblique  angle  at  the  base 
of  the  layered  soil  system  from  an  underlying  uniform  half  space. 

The  fundamental  equations  of  motion  are  presented  and  partially 
solved  in  Section  3.2.  The  complete  solution  for  the  special  case  of  a 
half  space  with  a  free  boundary  is  presented  in  Sections  3.3  and  3.4. 
Additionally,  an  exact  solution  for  a  single  uniform  layer  over  a 
viscoelastic  half  space  is  developed  in  Section  3.5.  This  solution  will 
be  later  used  to  verify  the  discretized  method,  as  developed  in  Sections 
3.6  and  3.7,  for  a  general  multi-layered  system. 

3.2  Governing  Equations 

The  motions  created  by  incident  plane  body  waves  will  in  general 
involve  three  components  of  displacements,  however,  these  components  do 
not  vary  in  the  horizontal  direction,  y,  perpendicular  to  the  direction 
of  wave  propagation.  Hence  the  problem  involves  only  the  space 
coordinates  x  and  z.  The  coordinates  of  x  and  z  are  defined  as  shown  in 
Fig.  3.1. 

For  the  special  case  of  harmonic  excitation  the  equations  of  motion 
for  an  isotropic  viscoelastic  medium  may  be  written: 
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(3. l.b) 
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3x2  3z2 


(Dilation) 


(Laplace  Operator) 


and  ux,  uy,  u2  are  the  displacements  in  the  x,  y,  z  directions, 
respectively,  p  is  the  mass  density  and  M*  and  G*  are  the  complex 
moduli  introduced  earlier.  Equations  (3.1. a)  and  (3. l.b)  are  coupled. 
They  govern  the  motion  in  the  vertical  x-z  plane  while  Eq.  (3.1.C) 
governs  the  motion  perpendicular  to  the  x-z  plane.  The  waves  described 
by  Eqs.  (3.1. a)  and  (3. l.b)  will  be  discussed  first. 

General  Solutions 

Using  a  method  based  on  the  work  of  Helmholtz,  see  Ewing  et  al. 
(1957),  Eqs.  (3.1. a)  and  (3. l.b)  can  be  solved  by  expressing  the 
displacements  in  terms  of  displacement  potentials  $  and  \|>  as  follows: 

u  .34-1* 


x  3x  3z 

u  -  |4  ♦  3* 
z  3z  3x 

where  4,  satisfy  the  wave  equations: 


(3. 2. a) 

(3.2. b) 
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V20  -  -jr  (3.3.b) 

V  3t 

■ 

The  potential  functions  4>  and  are  associated  with  P  and  SV  wave 

motions,  respectively.  The  two  types  of  waves  propagate  through  the 

medium  with  velocities  V  and  V  ,  respectively. 

p  s 

Since  the  motions  are  assumed  to  be  harmonic  at  the  frequency  ui, 
the  wave  potentials  must  also  be  harmonic,  and  can  be  written: 

♦  -  »  (3. 5. a) 

-  ¥  eiu)t  (3.5.b) 

where  ♦  and  ¥  are  complex  potential  amplitudes.  Substitute  of  Bqs. 
(3.S)  into  Bqs.  (3.3)  yields  the  time  independent  equations: 


V2*  +  k2  -  0  (3. 6. a) 

V2*  +  k2  -  0  (3. 6. a) 

3 

*  * 

where  k  »  u/V  and  k  ■  ai/V  are  complex  P-wave  and 
p  p  s  s 

S-wave  numbers,  respectively. 

Particular  solutions  to  Bqs.  (3.6)  corresponding  to  plane  waves 
propagating  in  the  positive  x-di recti on  with  the  complex  wave  number 
k  ■  kf  +  ik^  (k^  >  0,  kj  <  0)  may  be  found  by  separation  of 
variables: 

Let  *  -  f(i)  e'ilUt  (3.7.o) 

¥  -  g(z)  e"ikx  (3.7.b) 


Substitution  of  these  expressions  into  Bqs.  (3.6)  leads  to  the  following 
ordinary  differential,  equations: 

M  ♦  <k2  -  k2)  f  -  0  (3.8.a) 

dx2  p 
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^4  +  (k2  -  k2)  g  -  0 
dx^ 


(3.8.b) 


The  solutions  to  Eqs.  (3.8)  are: 

r  ik.z  -ik.z 

f(z)  *  A.  e  *  +  A_  e  . 


(3. 9. a) 


g(z)  -  B1  e  +  b2  e 


(3.9.b) 


where  k.  and  k,  are  the  solutions  with  positive  real  parts  to  the 

$  \p 

equations. 


k2  -  k2  -  k2 
4>  P  * 


(3. 10. a) 


2  2  2 

If  ■  If  *  |r 

>  s 


Introducing  the  notation 


k  •  k  sin  e  ■  k  sin  f 
P  * 


(3.10.b) 


(3.11) 


where  e  and  £  are  (not  necessarily  real)  angles  which,  as  will  be  shown 

in  connection  with  Bq.  (3.23),  are  related  to  the  direction  of  wave 

propagation,  the  wave  numbers  in  Bqs.  (3.10)  may  be  written 

k.  *  k..  cos  e  ■  k  cot  e  (3. 12. a) 

9  P 

k.  »  k  cos  f  »  k  cot  f  (3. 12. b) 


Substitution  of  Bqs.  (3.9)  and  (3.10)  into  Bqs.  (3.7)  yields  the 
following  solution  to  the  original  Helmholtz  equations,  Bqs.  (3.6): 


(ik.z  -ik.z  v  .. 

A^w)  e  *  +  A2«d)  e  *  j  e"lkx 

/  "lV\  -lkx 

f  -  (b1  (oj)  e  v  +  B2(tu)  e  w  )  e  lK* 


(3. 13. a) 


(3.13.b) 


A  further,  complete  expression  for  the  displacement  amplitudes  may  be 
found  below,  Bq.  (3.31).  The  surface  amplitudes  follow  from  Bq.  (3.33). 


*  * - 
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Equations  (3. 13. a)  and  (3.13.b)  are  obviously  of  the  same  form.  Hence, 
only  solutions  corresponding  to  the  first  equation  will  be  discussed  in 
full  detail.  These  solutions  correspond  to  P-waves. 

P- Waves 

The  potential  function  in  Eq.  (3. 13. a)  can  be  written  on  the  fora 

*  -  *  +  *  (3.14) 

where 

-i(K  *7)  — A  *r  -iP  »r 

♦  -  A  e  -  A  e  •  e  n  ;  n  -  1,2  (3.15) 

n  n  n 

In  this  notation  r  is  the  location  vector,  r  ■  xx  +  zz,  where  x  and  z 

are  unit  vectors  on  the  x-  and  z-axes,  respectively,  and  K  is  the 

n 

vector 


K  «  kx  +  (-l)Vz  -  P  -  iA  »  n  -  1,2  (3.16) 

n  9  n  n 

The  real  vectors  P  ■  Re(K  )  and  A  •  -lm(K  }  are  called  the 
n  n  n  n 

propagation  vector  and  the  attenuation  vector,  respectively. 

The  last  fora  of  Eq.  (3.15)  shows  that  $  represents  a  plane 

n 

wave  with  the  potential  aaplitude  A^  (at  the  origin)  which  propagates 

in  the  direction  P^  with  the  attenuation  factor  exp(-|An|)  per 

unit  length  in  the  direction  of  the  attenuation  vector  A  .  Both  of 

n 

these  vectors  lie  in  the  xz-plane.  A  wave  for  which  the  propagation 

vector  and  the  attenuation  vector  do  not  coincide  is  called  an 

inhonogeneous  wave.  For  such  plane  waves  the  amplitude  will  vary 

exponentially  along  the  wavefront,  see  Borcherdt  (1973) . 

The  angle  e  from  the  z-axis  to  the  propagation  vector,  F  may 
n  n 

be  determined  by  considering  the  vector  product 


A 

Z  X 


X 

0 

Re(k) 


y 

0 

0 


z 


1 

(-l)n  Re (k  ) 


Re (k) y 


(3.17) 
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where  y  is  a  unit  vector  on  the  y-axis.  This  vector  product  also 


satisfies 


z  x  P  ■  1*1  •  |P  I  sin  e  y  •  |P  I  sin  e  y 
n  1  1  1  n'  n  1  n'  n 


thus. 


sin  e  ■ 
n 


(3.18) 


(3.19) 


Jr«2 (k)  +  Re2 (k^) 


Similarly,  by  considering  the  "dot"  product  $*P 


(-1)"  Re(k , ) 


I  1  "  ■ "  ■ 

yj Re2(k)  +  Re2(k  ) 


(3.20) 


which,  since  Re(k.)  >  0  and  Re(k)  >  0,  shows  that  n  -  1  corresponds 

9 

to  a  wave  propagating  generally  upwards  while  n  ■  2  corresponds  to  a 
wave  propagating  generally  downwards,  see  Fig.  3.1. 

The  same  procedure  applied  to  the  vector  product  x 
yields  the  following  expression  for  the  angle,  a^,  from  the 
propagation  vector  to  ttv*  attenuation  vector 


sin  a  *  (-1) 
n 


Re (k , ) Im(k)  -  Im(k.)Re(k) 
_ 9  _ $ _ _ 

y  |ke2  (k)  +  Re2(k  )]  [lm2(k)  +  Im2(k.)] 


(3.21) 


This  expression  shows  that  a  homogeneous  wave  ( a ■  0)  occurs  if 


and  only  if 


Im(k) 


Im(k. ) 


Re (k)  Re (k. ) 

$ 


(3.22) 


But  this  implies,  by  Eqs.  (3.10)  that  arg(k)  ■  arg(k^)  ■  arg(k^), 
and  also  that  the  angle  e  in  Eg.  (3.11)  is  real  for  homogeneous  waves. 
Furthermore,  for  this  case  Eqs.  (3.19)  and  (3.20)  reduce  to 


sin  e„  »  sin  e 
n 

cos  e_  -  (-1) n  cos  e 
n 


(3. 23. a) 
(3.23. b) 
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Thu*,  •  is  siaply  the  angle,  «2,  which  the  downgoing  wave  fora*  with 
the  z-axis,  and  the  corresponding  angle  for  the  upgoing  wave  is 
•  w  -  e. 

SV-Wavc8 

A  siailar  vectorial  study  of  the  SV-wave  solutions  expressed  by  Bq. 

(3.13.b)  will  show  that  all  of  the  foraulas,  Eq.  (3.14)  to  Eq.  (3.23) 

are  also  valid  for  this  case  provided  k.,  k  ,  9  ,  A  and  e 

9  p  n  n  n 

are  replaced  by  k. ,  k f  ,  B„  and  f  ,  respectively, 
if  8  n  n  n 

SH-Wavea 

The  notions  described  by  Bq.  (3.1.C)  are  called  SH- waves.  They 
involve  displaoeaents  only  in  the  indirection.  The  solution  for 
harmonic  notions  propagating  with  the  wave  number  k  in  the  positive  x- 
direction  is 

uy  -  0y  eiut  (3. 26. a) 

where 


/ 

(Cl* 


(3.26.b) 


The  wave  number  k^  is  as  defined  by  Bq.  (3.10.b),  and  and  C2 
are  to  be  considered  arbitrary  complex  constants.  As  for  the  case  of 


the  P-  and  SV-wave s  discussed  above,  the  terns  of  this  solution  nay  be 


written  in  a  vector  notation  siailar  to  Bq.  (3.13)  and  (3.14).  However, 


homogeneous  incident  SH-waves  do  not  give  rise  to  inhoaogeneous 
reflected  waves,  and  the  simple  notation  of  Bq.  (3.2f)  will  suffice  for 
wave  fields  considered  in  this  dissertation.  The  wave  numbers  for 


homogeneous  SH-waves  satisfy  the  condition  stated  by  Bq.  (3.11),  l.e., 
k  ■  ka  sin  f,  where  f  is  the  angle  which  the  propagation  vector  foras 


with  the  s-axis 
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3.3  The  Viscoelastic  Half  Space 

The  following  study  of  the  solution  for  the  esse  of  a  homogeneous 
harmonic  body  wave  obliquely  incident  on  the  free  surface  of  a  uniform, 
viscoelastic  half  space  will  provide  an  insight  into  the  nature  of  the 
more  complicated  solutions  for  multi-layered  systems  to  be  discussed  in 
later  sections.  As  shown  in  Fig.  3.1,  an  obliquely  incident  txamogeneous 
SV-vave  will  in  general  result  in  two  reflected  waves  of  different  types 
(mode  conversion).  Furthermore,  in  certain  cases,  depending  on  the 
incident  angle  and  the  material  properties,  one  of  the  reflected  waves 
may  be  inhomogeneous,  Borcherdt  (1971)  and  Cooper  (1967). 

Incident  SV-Waves 

Only  the  specific  case  of  an  incident  SV-wave  will  be  discussed  in 
detail.  However,  the  method  presented  is  also  applicable  to  the  case 
of  an  incident  P-wave.  Since  no  incident  P-wave  exists  for  the  case 
studied,  the  form  of  the  solution  is  given  by  Eq.  (3.13)  with  A^  ■  0. 

The  incident  SV  wave  is  represented  by  the  term  exp(ikz). 

The  incident  angle,  f  ■  -fg,  shown  in  Fig.  3.1,  is  real  since 
the  wave  is  assumed  to  be  homogeneous.  The  relations  in  Eqs.  (3.23) 
immediately  implies  that  the  reflected  SV-wave,  represented  by  the  term 
B2  exp(ik^z)  is  homogeneous  and  forms  the  angle  f2  *  f  with  the 
z-axis.  Inversion  of  Eq.  (3.11)  and  multiplication  by  u  yields  the 
complex  form  of  Snell's  Law 

V*  »  V*/sin  e  -  V*/sin  f  (3.27) 

a  p  a 

* 

where  V  •  u/k  is  the  apparent  eomplex  phase  velocity  along  the  free 
8 

surface.  Hence, 

V_ 

-  -f  sin  f 

V 

s 


sin  e 


(3.28) 
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which  shows  that  the  reflected  P-wave,  corresponding  to  the  last  term 

in  Eg.  (3.13a) t  will  always  be  inhomogeneous  if  6#  *  8pi  and  that 

even  if  0  ■  $  a  homogeneous  reflected  P-wave  will  only  occur  when 

s  p 

f  is  smaller  than  a  certain  critical  angle,  f  ,  defined  by 

cr 


sin  f  ■  V*/V*  (3.29) 

cr  s  p 

since  sin  e  >  1,  for  f  >  f  . 

cr 

The  amplitudes  and  B2  of  the  reflected  P-  and  S-waves, 
respectively,  can  be  computed  from  the  condition  that  no  stresses  occur 
at  the  free  boundary.  Introducing  the  notation 

a  ■  cot  e;  b  ■  cot  f  (3.30) 

where  b  is  real  while  a  could  be  complex,  and  remembering  Eqs .  (3.12), 
the  displacement  amplitudes  are  by  Bgs.  (3.2),  (3.5),  (3.7),  and  (3.9) 


where 


(3.32) 


(3.33) 


The  strain  field  can  be  obtained  by  differentiation  of  Eg.  (3,31)  and 
the  stresses  then  follows  on  the  basis  of  Hooke's  Law.  This  leads  to 
the  following  expression  for  the  normal  stress,  a,  and  shear  stress,  t, 


at  z  •  0 
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*  1 

2a 

-(1-b2)  -2a  -(1-b2) 

U 

«  G  k2 

0 

(1-b2) 

2b  (1-b2)  -2b  _ 

fA, 


(3.34) 


Since  these  stresses  must  vanish,  A^  ■  0,  and  B^  is  known,  this 
immediately  leads  to  a  system  of  linear  equations  from  which  A^  and 


B2  nwy  be  determined. 


-2a 


- (1-b  ) 
-2a  _ 

The  solution  is: 


(1-b2) 


(3.35) 


.  .  -4b  (1-b) 

2  .  w  ^  „  .2, 2  "l 

4ab  +  (1-b  ) 

a  .  iab_z_ilzbif 

5  5  “l 
4ab  +  (1-b*)*  A 

This  solution  has  the  properties: 

A2  0  and  B2  ♦  -lj  for  f  0 
and 


(3.36) 


(3.37) 


A2  0  and  B2  **■  1}  for  f  45° 

Thus,  no  reflected  P-waves  occur  at  these  angles  of  incidence. 

Furthermore,  by  Bq.  (3.33),  no  vertical  displacement  occurs  at  the 

surface  for  f  ■  0;  and  no  horizontal  displacement  occurs  for  f  *  45°. 

Surface  displacement  amplitudes  for  other  angles  of  incidence  can  be 

computed  from  Bq.  (3.33). 

Solutions  for  the  elastic  case  have  been  published  by  Knopoff 

et  al.  (1957)  and  Meissner  (1965)  for  both  incidence  SV-  and  P-waves. 

Figs.  3.3  and  3.4  shows  ratios  between  surface  amplitudes  and  the 

horizontal  component  of  amplitude  of  a  vertical  incident  SV-waves.  The 

corresponding  particle  motions  are  shown  in  Fig.  3.2.  These  figures, 

which  are  actually  valid  for  a  damped  half  space  as  long  as  6  »  B  , 

*  P 


indicate  several  interesting  features: 
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3.4  Vertical  Component  of  Surface  Kotl 
Incident  SV  waves  -  Elastic  Ha If spec 


1.  For  angles  of  Incidence  less  than  the  critical  angle,  see 
Eg.  (3.29),  a  is  real  and  the  horizontal  and  vertical 


components  are  in  phase.  The  particle  notion  is  linear  and 
the  horizontal  component  of  motion  is  larger  than  the  vertical 
component  for  usual  values  of  Poisson's  ratio. 

2.  For  incident  angles  larger  than  the  critical  angle,  a  is 
complex. 

The  reflected  P-wave  is  inhomogeneous  and  the  reflected 
(homogeneous)  SV-wave  is  out  of  phase  with  the  incident  wave.  The 
horizontal  and  vertical  motions  are  90°  out  of  phase.  Thus,  the 
particle  motion  is  elliptical  (retrograde  for  f  <  45°  and  prograde 
for  f  >  45°) .  The  vertical  component  is  generally  larger  than  the 
horizontal  component. 

The  comparison  between  the  induced  vertical  and  horizontal 
amplitudes  of  surface  motion  for  an  incident  SV-wave  at  different 
angles  of  incidence  is  shown  in  Fig.  3.7.  The  dotted  lines  show  the 
cases  for  all  incident  angles  exceeding  the  critical  angle.  A 
singularity  is  found  at  45  degrees  because  the  horizontal  displacement 
is  0  for  all  values  of  Poisson's  ratio. 

Incident  P-Wave 

The  case  of  an  incident  P-wave  can  be  treated  by  the  technique 

used  above.  Snell's  law,  Eq.  (3.27)  is  also  valid  for  this  case,  but 
•  * 

since  V  is  always  larger  than  V  no  critical  angle  of  incidence 
P  s 

exists.  Surface  amplitudes  can  also  be  computed  from  Eq.  (3.33)  by 
setting  B^  *  0  and  assuming  that  A^  is  given.  Figures  3.5  and  3.6 
show  the  normalized  vertical  and  horizontal  amplitudes  of  P-wave  at 
different  angles  of  incidence  for  a  range  of  Poisson's  ratio.  Both 
displacement  amplitudes  are  normalized  by  the  amplitude  of  the 
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.  3.5  Vertical  Component  of  Motion  for  Inc 
P  waves  (after  Knopoff  at  al.,  1957) 


Ki«.  3.7  Displacement  Ratio  of  Surface  Fig.  3.5  Displacement  Ratio  of  Surface 

Motion  for  incident  SV  wave  Motion  for  Incldent  P_*ave 
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vertically  incident  P-wave.  The  features  of  the  notion  are  summarized 
as  follows: 

1.  For  a  medium  with  Poisson's  ratio  greater  than  0.2,  no 

significant  difference  is  seen  in  the  different  amplitude 

curves  and  an  approximation  through  a  function  U  /v 

z  o 

■  cos  e  is  sufficient  for  most  practical  purposes,  Meissner 

(1965).  V  denotes  the  vertical  amplitude  of  vertically 
o 

incident  P-waves. 

2.  The  induced  horizontal  components  start  from  0  at  vertical 
incidence,  and  increase  practically  linearly  to  35  degrees 
incidence.  The  maximum  horizontal  amplitude  is  about  52 
percent  of  Vq  at  52  degrees  incidence  for  a  medium  with  a 
Poisson's  ratio  of  0.4;  and  is  about  93  percent  of  Vq  at  65 
degrees  incidence  for  a  medium  with  Poisson's  ratio  of  0.2. 

The  ratio  of  the  horizontal  to  the  vertical  amplitude  of  surface 
motion  due  to  incident  P-wave  is  shown  in  Fig.  3.8.  The  angle  of 
incidence  has  a  strong  effect  on  the  displacement  ratio  for  materials 
with  a  low  Poisson's  ratio,  while  for  materials  with  a  high  Poisson's 
ratio  the  effect  is  relatively  small. 

3.4  Single  Layer  over  Half  Space 

In  order  to  verify  the  algorithm  of  the  computer  programs  SITE  and 
LOVE  which  will  be  described  in  Section  3.7,  the  boundary  value  problem 
for  a  viscoelastic  uniform  layer  overlying  a  viscoelastic  half  space  is 
solved  analytically.  The  amplification  of  oblique  SV-  or  P-waves 
incident  to  the  base  with  various  incident  angles  were  examined.  The 
effects  of  S-wave  velocity  ratios  between  the  layer  and  the  half  space 


on  response  were  also  studied 


3.4.1.  Solutions  to  Boundary  Value  Problem 


The  structural  model  shown  in  Fig.  3.9  consists  of  two  uniform 

isotropic  viscoelastic  media.  A  layer  with  the  thickness  H,  and  the 

properties  p* ,  G*'f  M*  ,  V*  ,  V*  ,  and  6^  overlies  a  half  space 

*  *  *  * 

with  the  properties  p,  G  ,  M  ,  Vgf  v^,  Bg,  and  0^. 

Let  ♦,  and  ♦' ,  4"  be  the  displacement  potentials  of  the  P-and 
SV-waves  in  the  bottom  half  space  and  stratified  layer.  The  general 
expressions  of  and  are  denoted  by 


.  ,,  iakz  .  -iakz,  -ikx 

4>  «  (A  e  +  A2  e  )  e 


ibkz  „  -ibkz,  -ikx 
Y  *  (Bx  e  +  b2  e  )  e 


where 


{'VV2  -  !}1/2  •  0  ■  {<V>»'2  -  J}: 


<3. 38. a) 


(3.38.b) 


(3.39) 


*• 

ir  kz 

-ir kz. 

-ikx 

-  (C  e 

+  D 

e 

) 

e 

iskz 

-iskz. 

-ikx 

4" 

*  (E  e 

+  F 

e 

) 

e 

where 


(3. 40. a) 

(3. 40. b) 

(3.40.C) 


t  •  {<V>*12  -  l}l/j  .  s  -  {|V*/V*’|2  -  lj1/2  (3.40. 

in  which  A^,  B2,  C,  D,  E,  F  are  arbitrary  constants  to 

be  determined  by  the  boundary  condition.  On  the  basis  of  assumptions 
similar  to  those  made  in  the  previous  section: 

1.  Only  one  incoming  wave  is  incident  to  the  base  boundary 


2.  The  incident  angle  is  real 

3.  V*  «  V*/sin  f 

a  s 

it  can  be  concluded  that  if  only  an  SV-wave  is  incident  to  the  base 
boundary,  the  coefficient  "b”  will  be  the  only  real  quantity  and  a,  r, 
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Fig.  3-9  Structural  Modal  -  Inclined  Body  Have 
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and  s  will  be  coaplex.  if  a  P-wave  is  incident,  then  b,  r,  and  a  will 
be  coaplex  in  general  and  the  coefficient  "a"  will  be  real. 

The  displacement  components  can  be  determined  from  Eg.  (3.2)  and 
stresses  can  also  be  obtained  from  Hooke's  Law  by  the  substitution  of 
appropriate  quantities  for  the  half  space  and  the  surfacial  layer.  The 


boundary  conditions  are: 
at  2  -  0:  Ux  -  Ux 


T  ■  T 
xz  xz 

at  z  ■  -H:  Txz  -  0  , 


U  ■  U 
z  z 


zz 


zz 


zz 


Incident  SV-Wave  From  the  Base 


By  setting  *  0  and  normalizing  all  coefficients  by  B^,  a 

set  of  linear  equations  may  be  formed  in  terms  of  the  system  geometry , 

system  properties,  and  normalized  coefficients.  The  linear  equation  is: 

□0{q  }  •  R  (3.41) 

s  s 


where 


(0S}- 

<A',  B', 

C',  D',  E' 

,  *’>T 

<-b,  1, 

2  * 

-(b  -1)C  , 

-2bG* ,  0,  0>T 

and  [p]  - 

1 

-b 

-1 

-1 

-s 

s 

-a 

-1 

-r 

r 

1 

1 

* 

-2aG 

2  * 

(b  -1)G 

*  < 

-2rG 

•  • 

2rG 

2  *' 
-(s  -1)G 

2  *' 
-(s  -1)G 

2  * 

- (b  -1)G 

* 

-2bG 

(s2-1)G*' 

*  • 

-2sG 

*• 

2sG 

0 

0 

-  ihkr 
-2re 

(s2-l)e’ik8h 

0 

0 

-(s2-l)e-ikhr  -(a2-l)eikhr 

2se’iksh 

-2Mik,h 

(3.42) 

•  I 

From  this  equation  the  normalized  coefficients  Bj,  C',  D', 

£' ,  F'  can  be  found  in  terms  of  given  parameters  in  CpD  and  {r#}. 
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Incident  P-Wave  Proa  the  Base 

Again,  a  real  value  of  an  angle  of  incidence  may  be  specified, 
then,  by  setting  B^  -  0,  and  normalizing  all  coefficients  by  A^s 

*  {*P}  <3-43> 

in* which 

|Qp(  «  <A^,  B',  C’,  O',  E',  P'>T 

and 

jR  )  -  <-l,  -a,  -2aG*,  (b2-l)G*,  0,  0>T 
Displacements  and  Stresses 

For  any  wave  field  expansion  of  Eq.  (3.2)  gives  the  normalized 

displacements  for  the  surface  layer: 
u.  «  y  ,ei  (oit-kx) 

(3.44) 

w'  =  Wei(u)t‘kx) 

where 


U’ 

«  -ik 

|(C' 

ikrz  ... 
e  +  D' 

-ikrz 

e 

W' 

-  ik 

|r  (C' 

ikrz 

e  -  D' 

-ikrz 

e 

.  iksz  .  iksz 

)  -  (E'  e  +  F'  e 


>) 


(3.45) 


and  the  stress  components 
i  (wt-kx) 


T'  *  t  e 
xz  - 


o'  ■  o  e 
xz  - 


i  (wt-kx) 


(3.46) 


where 


#* 


.  '  .<i,  .  irkz  „ .  -irxz,  ,  i  ibkz 

2  ■  G  k  {2r  (C •  e  -  D'  e  )  +  (s  -1) (E'  e  ♦  F'  e 


-irkz. 


iskz 


-iskz 


)} 


(3.47) 


#• 


£  •  G  k2{- (s2-l) <C •  eirkz  +  K'  e_lrkz)  -  2s(E'  eiikz  -  F'  ei8kz) } 

The  displacements  and  stresses  at  the  base  boundary  due  to  a  specified 
incident  wave  take  the  same  forms  as  shown  in  Egs.  (3.33)  and  (3.34) 
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3.4.2  Numerical  Example! 

A  aodal  consisting  of  a  siaple  uniform  layer  overlying  a  half  space 
was  used  to  investigate  the  characteristics  of  displacements  at  the  free 
surface  and  the  base  boundary  for  the  case  of  SV-  or  P-waves  incident 
at  different  angles.  The  study  also  included  a  study  of  the  effects  of 
the  incident  angle  on  amplification  of  both  steady  state  and  transient 
motions.  The  influence  of  the  variation  of  shear  wave  velocity, 

Poisson's  ratio,  and  damping  ratio  between  surface  layer  and  half  space 
were  also  studied. 

The  basic  system  is  shown  in  Fig.  3.10.  The  parameters  had  the 
values 

1 


Surface  Layer 

Half  Space 

Thickness  (ft) 

128 

infinite 

Unit  weight  (pcf) 

125 

162.5 

Shear  wave  velocity  (fps) 

1000 

1200  -  8000 

Poisson's  ratio 

0.1  -  0.45 

0.1  -  0.45 

Damping  ratio 

0.05  -  0.12 

0.02  -  0.05 

The  fundamental  fixed  base  undamped  frequencies  of  the  layer  were 


1.95  Hz  and  3.38  Hz  for  S-  and  P-waves,  respectively.  Computations 
were  performed  at  a  frequency  of  1  Hz.  Thus,  the  shear  wave  length  in 
the  surface  layer  was  1,000  ft.  The  ratio  of  wave  length  to  the 
thickness  of  the  layer  was  about  7.8  and  the  ratio  of  densities  between 
half  space  and  the  surface  layer  was  1.3. 

Displacements  at  Laver  Surface  and  Laver  Base 

The  horizontal  displacements  and  vertical  displacements  at  the  free 
surface  due  to  a  harmonic  inclined  SV-wave  was  computed  for  an  elastic 
system  with  Poisson's  ratio  0.25.  The  results  are  shewn  in  Pig.  3.10. 

In  order  to  see  the  influence  of  the  incident  angle,  both  components  are 
normalized  by  the  displacement  at  ground  surface  due  to  a  vertically 
propagating  shear  wave.  The  results  are  very  similar  to  those  shown  in 
Section  3.3  for  a  viscoelastic  half  space.  The  horizontal  component  has 
a  sharp  peak  at  the  critical  angle  whereas  the  vertical  component  at 
this  point  forms  a  downward  sharp  cusp.  For  most  realistic  choices  of 
system  parameters,  the  horizontal  amplitudes  of  inclined  SV-waves  are 
maaller  than  those  of  vertically  propagating  shear  waves  except  at  a 
very  narrow  range  of  angles  of  incidence  around  the  critical  angle. 

The  vertical  amplitudes  of  SV-waves  are  zero  at  normal  incidence  and 
linearly  increase  to  40  percent  of  the  amplitudes  of  vertically 
propagating  shear  waves  at  30  degrees  of  incidence.  The  effects  of  the 
shear  wave  velocity  of  the  half  space  on  the  displacements  are 
insignificant,  as  can  be  seen  from  Fig.  3.10  which  contains  results  for 
four  different  values  of  the  velocity  ratio. 

Figure  3.11  shows  similar  results  for  the  case  of  incident  P-waves. 
Both  components  are  normalized  to  the  ground  surface  amplitude  of  a 


► 


Amplitude  of  vortical  P-wave 
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vertically  incident  P-wave.  All  the  vertical  components  of  the 
inclined  P-wave  are  mailer  than  those  corresponding  to  vertical 
incidence.  No  difference  can  be  seen  in  the  amplitude  curves  for  the 
four  values  of  shear  wave  velocity  ratios.  This  effect  is  probably  due 
to  the  relatively  long  wave  length  as  ccaipared  to  the  thickness  of  the 
surfacial  layer.  The  horizontal  components  of  the  inclined  P-wave  are 
more  sensitive  to  the  shear  wave  velocity  of  the  half  space.  The 
greater  the  stiffness  of  the  elastic  half  space,  the  larger  the 
horizonal  displacements  due  to  inclined  P-waves.  The  greatest 
differences  in  horizontal  amplitudes  occur  for  angles  of  incidence 
between  50  and  70  degrees. 

Figure  3.12  illustrates  similar  studies  of  normalized 
displacements  at  the  top  of  the  half  space  for  the  case  of  a  shear  wave 
velocity  ratio  equal  to  4.  Both  results  indicate  that  if  the  angle  of 
incidence  is  less  than  30  degrees,  which  should  be  the  practical  case; 
the  vertically  propagating  body  waves  produce  larger  motions  than  the 
inclined  body  waves. 

Amplification  at  Ground  Surface 

A  study  was  made  of  the  variation  of  amplification  of  a  harmonic 
SV-wave  arriving  at  a  wide  range  of  incident  angles.  The  frequency  for 
this  harmonic  SV-wave  was  arbitrarily  chosen  as  1  Hz.  Site 
amplification  at  ground  surface  versus  angles  of  incidence  have  been 
plotted  for  four  different  values  of  the  S-wave  velocity  ratio  as  shown 
in  Fig.  3.13.  For  angles  of  incidence  less  than  the  critical  angle, 
the  horizontal  components  at  the  ground  surface  are  about  1.4  times  the 
horizontal  amplitudes  of  SV-wave  at  the  base  for  the  cases  of  an  S-wave 
velocity  ratio  greater  than  2.  The  vertical  components  at  the  normal 
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ANCLE  OF  INCIDENCE  -  Deg.  ANCLE  OF  INCIDENCE  -  Deg. 

Fig.  3.12  Horizontal  and  Vertical  Component  of  Motion  at  Top  of  Bedrock  for  Incident  Body  naves 


Amplification  at  Ground  Surface  for  Incident  SV  waves 


incidence  ace  zero  and  increase  linearly  to  about  70  percent  of  the 
horizontal  amplitude  of  the  vertical  shear  wave  at  30  degrees  of 
incidence.  At  the  critical  angle,  vertical  components  show  a  very 
sharp  trough  because  an  unusual  amplitude  spike  appears  at  this  angle 
as  seen  in  Pig.  3.12.  The  amplification  at  seme  particular  angle  of 
incidence,  say  43  degrees  in  this  case,  can  not  be  defined  since  the 
horizontal  amplitude  at  the  base  approaches  zero.  Fortunately,  this 
case  will  not  appear  in  practice  because  the  angles  of  incidence  are 
normally  less  than  the  critical  angle,  and  also  this  phenomenon  does 
not  occur  in  damped  systems.  The  amplifications  for  inclined  P~waves 
incident  to  the  base  are  shown  in  Pig.  3.14.  Both  components  are 
normalized  by  the  vertical  amplitudes  of  the  base  at  the  corresponding 
angles  of  incidence.  For  this  case  the  shear  wave  velocity  ratio  is 
not  a  significant  factor  if  the  angle  of  incidence  is  less  than  the 
critical  angle. 

The  variation  of  amplification  with  frequency  for  an  undamped 
system  is  shown  in  Pig.  3.15  and  for  angles  of  incidence  of  0,  10,  20, 
and  30  degrees.  The  case  of  normal  incidence,  shown  as  a  solid  curve, 
is  the  well  known  case  of  a  vertically  propagating  shear  wave.  Only 
horizontal  amplitudes  occur  in  this  case.  The  peak  amplifications 
occur  at  1.95  Hz,  5.86  Hz,  9.77  Hz,  13.67  Hz  etc.,  which  are  the 
natural  frequencies  of  the  system.  For  inclined  waves,  coupling 
effects  occur,  and  vertical  components  are  induced.  Three  other  cases 
were  shown  in  the  same  plot.  Some  effects  of  coupled  P  wave  motion  may 
be  seen  at  frequencies  near  3.4  Hz  and  10.2  Hz,  which  are  the  first  and 
second  natural  frequency  of  P-waves  of  the  site. 


ANGLE  OF  INCIDENCE  -  Deg.  ANCLE  OF  INCIDENCE  -  Deg. 

Fig.  3.14  Anpllf lcatlon  at  Ground  Surface  for  Incident  P  waves 
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Computations  were  performed  for  the  same  system  assuming  a  2 
percent  damping  in  the  half  space,  and  5  percent  damping  in  the  layer. 
The  results  are  shown  in  Fig.  3.16.  The  resonant  peaks  for  all  modes 
are  damped  out  significantly  compared  to  the  undamped  cases.  However, 
while  the  effect  of  incident  angles  on  the  response  is  not  important  on 
the  horizontal  component,  it  is  quite  significant  on  vertical  component 
of  SV-waves  as  was  found  in  the  undamped  case. 

Effect  of  Damping  Contrast 

As  shown  above  the  effect  of  uniform  damping  is  to  reduce  the 
surface  response  significantly.  In  order  to  study  the  effect  of  a 
contrast  in  damping  ratio  between  the  surface  layer  and  the  underlying 
half-space  surface  amplification  factors  were  computed  for  the  special 
case  of  an  SV-wave  with  an  incident  angle  of  10  degrees  for  the 
following  three  choices  of  damping  ratios. 

Surface  Layer  Half  Space 
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The  result  of  these  computations  are  shown  in  Fig.  3.17.  As 
expected  the  horizontal  amplification  functions  were  found  to  be 
similar  to  those  computed  for  vertically  propagating  shear  waves, 
compare  Fig.  3.16.  However,  some  effects  were  observed  in  the 
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amplifications  of  vertical  notions.  In  particular,  it  is  interesting 

to  observe  that  Systems  3  (which  6g  ■  12%,  8p  ■  6%  in  the  layer) 

produced  nearly  as  big  vertical  notions  as  Systen  1  (with  8  -6 

s  p 

■  5%  in  the  layer).  Thus,  it  appears  that  horizontal  amplification  is 
governed  by  8_  while  vertical  amplification  is  governed  by  8  . 

8  p 

This  observation  is  perhaps  not  too  surprising  if  one  observed  that  for 
Systen  3  the  peaks  occur  at  the  natural  frequencies  for  vertically 
propagating  P- waves  in  the  site. 

Effect  of  Poisson’s  Ratio 

Studies  were  also  performed  on  the  above  undamped  two-layer 
systems  to  determine  the  influence  of  Poisson's  ratio.  Again,  a  SV 
wave  with  an  incident  angle  of  10  degrees  was  assumed. 

Figure  3.18  shows  the  effect  of  varying  Poisson's  ratio  in  the 
half-space  (it  was  set  to  0.25  in  the  layer)  and  Pig.  3.19  shows  the 
effect  of  varying  Poisson's  ratio  on  the  surface  layer  (it  was  set  to 
0.25  in  the  half-space).  As  expected  variations  in  Poisson's  ratio  has 
only  small  effects  an  the  horizontal  response  while  it  significantly 
influenced  the  vertical  response.  This  of  course  is  due  to  the  strong 
effect  of  Poisson's  ratio  on  the  P-wave  velocity. 

3.5  Multi-lavered  Half-space  for  SV  and  P  Waves 

The  propagation  of  body  waves  in  multi-layered  systems  is  of 
fundamental  interest  in  seismology  and  has  been  studied  by  Thomson 
(1950),  Matsumoto  (1953),  Haskell  (1960,  1962),  Phinney  (1964),  Hannon 
(1964),  Teng  (1967),  and  Bakun  (1970).  Their  methods  of  analysis  were 
essentially  based  on  the  Thomson-Haskell's  matrix  formulation.  The 
formulation  uses  the  idea  of  displacement  potential  theory  and  only 
considers  harmonic  wave  propagation  in  an  undamped  elastic  media. 
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Silva  (1976)  extended  Thomson-Haskell's  formulation  to  include  damping 
in  a  layered  system  overlying  an  undamped  elastic  half  space. 

The  above  theories  lead  to  a  very  complicated  equation  of  motion 
in  which  the  wave  number  k,  enters  through  transcendental  functions  in 
a  form  similar  to  Eq.  3.41.  As  shown  by  Lysmer  (1969)  and  Waas  (1972) 
a  much  simpler  equation  of  motion,  see  Eq.  (3.50),  can  be  obtained  by  a 
discretized  representation  of  the  layered  half  space.  In  this  method 
it  is  assumed  that  displacements  vary  linearly  between  layer 
interfaces.  As  will  be  discussed  this  assumption  imposes  certain 
restrictions  on  how  thick  individual  sublayers  can  be  chosen.  A  major 
advantage  of  the  method  is  that  it  leads  to  a  uniform  treatment  of 
inclined  body  waves  and  surface  waves  and  it  will  be  used  extensively 
below.  In  the  original  applications  of  the  method,  Lysmer  (1969)  and 
Haas  (1972) ,  the  underlying  half  space  was  considered  to  be  rigid. 

This  lead  to  a  method  for  analyzing  surface  waves  in  layered  systems  as 
will  be  discussed  in  Chapter  5. 

Udaka  (1975)  used  the  same  method  with  specified  motions  at  the 
half-space  surface  to  simulate  the  effect  of  traveling  waves.  In  this 
method  a  control  motion  at  the  half-space  surface  was  assumed  to 
propagate  horizontally  with  a  given  constant  phase  velocity.  Only  one 
component  of  motion  was  allowed  at  the  control  point.  The  method  does 
not  consider  interaction  between  the  layered  soil  systems  and  the 
underlying  half  space. 

This  interaction  is  properly  considered  in  the  present  work  on  the 
effects  of  inclined  body  waves  arriving  through  a  viscoelastic 
half  space.  The  study  will  show  that  except  for  the  case  of  normal 
incidence,  there  must  be  two  components  of  motion  at  each  layer 
interface.  The  following  cases  will  be  considered: 
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1.  SV- waves  at  oblique  incidence 

2.  P-waves  at  oblique  incidence 

3.  Pairs  o£  SV-and  P-waves  at  oblique  incidence 
3.5.1  Discretized  Formulation  for  Layered  System 

The  layered  system  considered  is  shown  in  Pig.  3.20.  It  consists 
of  n  homogenous/  isotropic  layers  over  a  half  space.  Ml  material 
properties  may  be  undamped  elastic  or  damped  viscoelastic. 

Assuming  linear  variation  of  displacements  within  each  layer  the 
displacement  of  any  point  in  the  system  for  a  plane  harmonic  wave 
travelling  in  the  horizontal  x-direction  at  a  given  frequency,  u>,  can  be 


written: 

6x  *  ««<*>  el(ut"kx)  (3. 48. a) 

«  -  ictfl(z)  ei(Wt-kx)  (3-48. b) 

•  z 

in  which,  for  z^  <  z  <  zJ+1 

Ux(2)  -  (zj+1  -  zJ/hj  U2j_1  +  (z  -  z^/hj  U2j+X  (3. 49. a) 

°z(s)  -  (Zj+i  “  z)/hj  U2j  +  (z  -  z^/hj  02j+2  (3.49. b) 


and  a  is  a  mode  participation  factor  that  can  be  found  from  the  given 
control  motion.  The  displacement  functions  Ux(z)  and  U^fz)  are 
interconnected  and  can  be  normalized  in  any  manner.  The  wave  number  K 
may  be  complex  expressing  both  the  phase  velocity,  V  ■  u/Re(k),  and 

A 

the  attenuation  factor,  exp  (-Zm(k)x). 

As  shown  by  Udaka  (1975)  the  equations  of  motion  for  the 
discretized  layered  system  is 

(DO  -  w2C(fl)  {u}  -  |  (3. 50. a) 

l  b/ 

whtrt 

DO  -  A  k2  ♦  C5]  k  +  [G]  (3. 50. b) 

In  these  equations  all  matrices  are  of  the  order  (2n  +  2)  x  (2n  +  2) 
and  the  last  two  terms,  P.  ,  of  the  load  vector  are  forces  at  the 

D 


interface  between  the  layered  system  and  the  halt  apace.  The  vector  0 


contains  (2n  *  2)  complex  displacement  amplitudes,  IT,  j  »  l,  2, 


(2n  +  2)  for  the  (n+1)  layer  interfaces  each  having  two  degrees  of 


freedom  U  and  U  . 

X  2 


The  banded,  symmetric  matrices  [A] ,  [§],  [G],  and  [m],  are  assembled 
from  sublayer  matrices  as  shown  in  Fig.  3.21.  The  submatrices,  shown  in 


Eqs.  (3.51)  to  (3.55),  are  formed  using  complex  shear  moduli  G  ^ , 

* 

complex  constrained  moduli  M.f  and  layer  thicknesses  h.  as  follows: 
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Figure  3.20  Kodel  of  Plane  SV  wave  Incidence 


Fig.  3.21  Structure  of  Matrices  [A],  [9],  [c],  and  [m] 
(For  the  case  of  SV  and  P  wave  of  oblique  incidence) 


The  layer  mass  matrix  M  can  be  chosen  as  consistent  mass,  or 
lumped  mass,  or  a  combination  of  both.  For  consistent  mass  formulation: 


tMc]J 


hj  1/6 


(3.54) 


For  lumped  mass: 


cVj  ‘  pj  hj 


(3.55) 


For  combination  of  consistent  and  lumped  mass: 


[ML  -  “cMi  * 11  -  vN: 


(3.56) 


where  a  is  a  fraction  between  0  and  1.  It  has  been  found  that  with 
c 

the  finite  element  method,  ■  0.5  gives  good  results  compared  with 

analytical  closed  form  solutions. 

The  same  investigation  also  showed  that  when  the  average  mass 

matrix  is  used  the  layer  thicknesses,  h^,  can  be  chosen  as  large  as 

h.  *  X  /5  ■  2flV  /(5w),  where  X  is  the  wavelength  of  shear 
j  s  s  s 

waves,  without  impairing  the  accuracy  at  the  discretized  method.  This 
compares  with  a  maximum  layer  thickness  of  Xg/8  when  either  the 
lumped  or  consistent  mass  formulations  is  used.. 

When  the  wave  number  k  becomes  zero,  i.e.,  the  apparent  wave 


velocity  approaches  infinity,  the  equation  of  motion  becomes 


|cg3  -  «*[«:}  ju}  .{°J 


(3.57) 


which  is  the  case  of  vertical  incidence  of  body  waves.  In  this  case 
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the  shear  wave  and  P-wave  are  completely  uncoupled  and  Eq.  (3.57)  can 
be  separated  into  two  simpler  matrix  equations;  one  for  each  wave  type. 

For  the  purpose  of  developing  the  base  boundary  equation,  it  is 
convenient  to  express  Eq.  (3.50a)  in  partitioned  form: 


in  which  the  suffix  "i"  denotes  quantities  for  the  layered  system,  the 
suffix  "b"  refers  to  the  base  boundary,  and  the  suffix  *c"  to 
interaction  between  the  layered  system  and  the  half  space.  For  a  given 
wave  number  k  and  frequency  u,  this  set  of  linear  equations  can  be 
further'  simplified  by  the  notation  DC  »CeJ  -  u  ft*]  to: 


“ _  _  - 
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or  [KJ  i0t)  .  [Kj  fo6}  .  (o) 
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Equation  (3.60)  yields 


(3.59) 

(3.60) 

(3.61) 


"  "Pd"1  Pc]  Kd  1  "  1#  2'  . '  2n  (3.62) 

and  substitution  of  this  equation  into  Eq.  (3.61)  gives  the  relationship 

between  the  displacements,  {u^},  and  the  forces,  {p^},  at  the  interface 
between  the  layered  system  and  the  half  space  (the  base  boundary) 


M  {%)  -  {Pb} 

where 


H  ■  -W1  W*1  [*J  *  W 


•» 

L 
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11 

L12 

21 

L22 

(3.63) 


It  should  be  noted  that  for  a  system  with  n  layers,  is  a  2n  by  2n 
matrix,  K.  is  a  2  by  2  matrix,  and  K  is  a  2n  by  2  matrix. 

D  C 
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3.5.2  Half  Space 

If  the  displacements  and  forces  at  a  specified  location  on  the 
base  boundary  induced  by  the  incident  wave  from  a  homogeneous 
viscoelastic  half  space  can  be  determined,  then  the  complex 
displacement  functions  {u^}  at  the  interfaces  between  each  layer  in 
the  system  can  also  be  easily  determined.  It  is  thus  necessary  now  to 
find  the  displacements  and  forces  at  the  base  boundary  due  to  a 
specified  incident  plane  wave  striking  at  the  base.  It  may  be  shown 
that  the  boundary  displacements  and  forces  take  the  same  forms  as 
Egs.  (3.33)  and  (3.34);  after  rearrangement,  they  are: 

lB 


{ub}  «  ik  [ot] 


and 


♦  G*  k2  [6'] 


(3.65) 


where 

1  -b" 

_  a  -1. 


H  • 

‘  2a 

(b2-lf 

.  - 

'-2a 

(b2-ir 

-  (b2-l) 

2b  . 

(b2-l) 

-2b  . 

and  Aj,  B1#  Aj,  and  are  complex  coefficients  representing 
wave  amplitudes  of  displacement  potentials. 

Elimination  of  boundary  forces  and  displacements  by  substitution  of 
Eqs.  (3.64)  and  (3.65)  into  Eq.  (3.63),  produces  two  linear  equations 
containing  four  unknown  coefficients.  As  discussed  earlier,  some 
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restrictions  and  assumptions  aust  be  aade  in  order  to  determine  the 
coefficients  and  to  obtain  normalized  boundary  displacements  and 
boundary  forces  at  the  base  boundary.  By  specifying  only  one  type  of 
body  wave  incident  to  the  base  by  a  given  real  incident  angle  and 
normalizing  two  reflection  coefficients  by  this  incident  coefficient, 
the  two  unknown  normalized  reflection  coefficients  can  be  determined 
from  the  following  equations: 


where 


cor1  • 

kl 

(3.66) 
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(D13  -  D14J 

(3.67) 
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in  which  D^,  d12,  d13,  D14,  D21,  D22,  D23,  and  D24  are  defined  as: 
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where  parameters  a  and  b  are  complex  parameters  as  defined  earlier  and 

1^,  l^2' L21'  L22'  ar*  elements  of  the  layer  stiffness  matrix 
defined  in  Eq.  (3.63). 

It  should  also  be  noted  that  the  coefficients  A2and  B2  are 
"Normalised  reflection  coefficients"  which  are  different  from  the 
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coefficients  *2  *nd  ®2  in  ®9S*  and  (3.65).  They  are  related 

to  the  reflection  coefficients  for  the  case  of  unit  incident 
ooefficeint.  It  should  also  be  noted  that  the  quantities  and  Q2 
will  be  redefined  from  case  to  case  for  different  incident  waves. 
Only  SV  Wave  Incidence 

Assuming  that  the  angle  of  incidence,  f,  is  real  and  that  the 

angle,  between  the  propagation  vector  and  the  attenuation  vector 

is  equal  to  zero,  one  may  obtain  k  ■  k  sin  f  as  shown  in  Eq.  (3.11). 

s 

Thus,  for  a  system  having  a  common  phase  velocity  and  attenuation 

* 

factor  along  the  boundary,  the  complex  wave  velocity,  V  is 

A 


*  *  *  ^V2 

V  *  u)/(k  sin  f)  *  V  /sin  f  *  V  e 

AS  S  A 


(3.68) 


in  which 
tan  £ 


28. Vl  -  62/(l  -  282)  28  for  small  8 

S  S  S  3 


a  s 

This  assumption  will  force  the  complex  parameter  ”b”  to  become  real  as 
long  as  Va  >  V_.  If  6D  *  8a  the  parameter  "a"  is  complex  and 


-  -  {(<//) 2  •  ei<V6a>  -  l}1/2 


where 


tan  6  -  28  Vl  -  82/(l  -  282)  28  for  small  8 

p  p  1  P  P  P 

If  6  ■  <5  *6,  the  solution  is  similar  to  the  elastic  case,  i.e., 

s  p 

both  "a"  and  "b"  are  real  but  the  wave  attenuates  in  the  direction  of 
propagation.  For  SV  wave  incidence,  A^  is  set  to  0  and  A2, 
can  be  determined  by  defining  and  Q 2  in  Eq.  (3.66)  by 


-<DU  +  D14> 

“(D23  +  °24> 

(3.69a) 

Once  two  normalized  coefficients  Aj,  B2  are  found  the  normalised 
displacements  and  boundary  forces  can  be  determined. 


Only  P-Wave  Incidence 

In  this  ease  is  set  to  zero  and  the  reflection  coefficients 

are  normalized  by  A^.  If  a  real  incident  angle,  e,  is  specified  we 
*  *  * 

have  V  *  V  /sin  e  and  k  *  u>/V  .  From  the  complex  Snell's  law, 
as  a 

*  *  * 

V  ■  v  /sin  f  ■  /V  /sin  e,  the  reflected  angle  f  can  be 

as  p 

determined.  Generally,  if  damping  due  to  the  shear  and  P  waves  are  not 
equal,  the  reflection  angle  of  the  SV  wave  will  be  complex  which  is 
different  from  the  elastic  case.  For  a  viscoelastic  half-space,  a 
phase  shift  at  the  boundary  base  will  generally  be  expected  regardless 
of  the  type  of  wave  incidence.  The  reflection  coefficients  can  be 
found  by  using  Eq.  (3.66)  by  replacing  and  Q2  as  follows: 


-(DU  +  °12> 

-<D21  +  °22) 

(3.69.b) 

SV-and  P-Waves  Obliquely  Incident 

It  is  assumed  that  the  amplitude  of  the  incident  P  wave  A^  is 
given  and  that  the  ratio,  n  *  o£  the  displacement 

coefficients  for  the  incident  SV  wave  and  the  incident  P  wave  is  also 
known.  With  these  assumptions  and  Q2  in  Eq.  (3.66)  will  be 
defined  as: 

°1  ■  -|BU  *  V  -  "<B1J  +  B14> 

Oj  ■  -102l  *  0j2)  -  n(D23  <•  D2<>  (3.69.C) 

Thus,  the  normalized  reflection  coefficients  can  be  determined.  In 
this  case,  the  reflection  angles  can  be  either  real  or  complex  and  will 
depend  on  the  damping  ratio  assumed  for  the  half-space  and  the 
specified  angle  of  incidence. 
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3.6  SH-Waves  in  Multi-layered  Half  Space 

Harmonic  notion  under  plane  atrain  conditions  in  a  semi-infinite 
layered  systen  due  to  an  inclined  SH-wave  incident  to  the  base  is  here 
considered.  The  model  for  plane  SV-and  P-wave  is  used,  and  is 
replotted  for  the  case  of  incident  SH-waves  as  shown  in  Fig.  3.22.  All 
displacements  are  perpendicular  to  the  x-z  plane  and  are  described  by 
i,..!,!.).1"**1  O.M) 

in  which  u>  is  the  circular  frequency,  and  a  is  the  mode  participation 
factor . 

The  equation  of  motion  take  the  same  form  as  Eq.  (3.50a)  but  the 
stiffness  matrix  00  Is  now  defined  as 

00  *  M  k2  +  [G]  (3.71) 


where  [A] ,  and  CG3  are  n  by  n  tridiagonal,  symmetric  matrices,  which 
consist  of  the  contributions  from  individual  layers  and  which  can  be 
conveniently  assembled  from  the  layer  submatrices  as  shorn  in  Fig.  3.23. 
The  submatrices  for  the  layer  j  are: 


P0 


j 


(3.72) 


(3.73) 


The  mass  matrix  (see  Eq.  3.56)  may  be  a  combination  of 


pL/3 

1/6] 

(consistent) 

Ni  *  P3  ’*3  , 

[l/6 

J 

(3.74) 

and 

fl/2 

l/6l 

(lumped) 

[Kt])  ■»)  hj 

. 

l/2j 

(3.75) 

Following  the  procedures  described  in  Section  3.5,  through  Eqs.  (3.58) 
to  (3.63),  the  relationships  between  the  boundary  force  and 


displacement  may  be  established.  On  the  other  hand,  the  boundary  force 
and  displaoeaent  induced  by  the  incident  SB-wave  frca  the  underlying 
half  space  can  be  derived  from  Eq.  (3.26)  in  the  way  similar  to  the 
case  of  SV- waves.  The  boundary  force  and  the  displaoeaent  are 
noraalized  by  the  amplitude  of  the  incident  wave  and  can  be  expressed 
in  teras  of  wave  number,  shear  modulus  and  and  unknown  reflection 
coefficient.  By  elimination  of  the  displacement  and  the  boundary  force 
between  the  layered  system  and  the  half  space,  the  normalized 
reflection  coefficient  can  be  found  and  the  boundary  force  and 
displacement  can  be  determined.  Accordingly,  the  amplitudes  of  each 
layer  can  be  easily  found  from  a  set  of  simple  linear  equations  for  any 
given  wave  number  and  frequency. 

3. 7  The  Computer  Programs  SITE  and  LOVE 

The  above  discretized  procedures  were  implemented  in  two  computer 
codes,  SITE  and  LOVE.  The  first  program  handles  the  cases  of  inclined 
P-  and  SV-waves.  It  also  handles  the  case  of  Rayleigh  waves  which  will 
be  discussed  in  Chapter  4.  The  second  program  handles  the  cases  of 
inclined  SH-waves  and  Love  waves  (also  to  be  discussed  in  Chapter  4) . 

Both  computer  codes  were  verified  against  the  exact  solutions 
provided  in  Fig.  3.5,  3.6,  3.10  to  3.12  and  3.15  to  3.17  using  a  model 
consisting  of  18  sublayers  to  represent  the  surface  layer.  The 
computed  data  points  agreed  with  the  exact  solutions  to  within  3 
significant  digits  and  have  not  been  plotted  in  the  above  figures. 

The  two  codes  can  handle  not  only  harmonic  motion  but  also 
transient  motions  by  the  Fourier  techniques  described  in  Chapter  5. 

3.8  Sunmary 


The  fundamental  theories  of  harmonic  inclined  body  waves 
propagating  in  a  viscoelastic  half  spaoe  and  a  layered  half  space  are 


82 


presented  in  this  chapter.  Baaed  on  analytical  solutions  developed  for 
a  nodal  of  a  single  layer  over  a  half  space,  several  miner ical  examples 
were  given  to  show  the  characteristic  of  amplitudes  at  the  free  surface 
and  the  interface  boundary  for  the  case  of  SV-  and  P-wave  incident  at 
different  angle.  The  study  of  the  effect  of  Incident  angle  on  the 
amplication  of  both  steady  state  and  transient  motion  was  also 
included.  The  effect  of  damping  and  Poisson's  ratio  contrast  (between 
the  layer  and  the  half  space)  on  site  amplifications  was  also 
investigated.  Numerical  examples  for  site  response  to  inclined  body 
waves  for  a  more  complicated  layered  half-space  system  will  be  given  in 
Chapter  5. 


CHAPTER  4 


PLANE  SURFACE  WAVES 


4.1  Introduction 

The  theory  of  plane  Rayleigh  and  Love  waves  propagating  in  an 
undamped  elastic  or  damped  viscoelastic  half  space,  a  layered  rigid 
base  system,  or  a  layered  half  space  has  been  well  developed. 

Basically,  the  available  theories  can  be  classified  into  continuum 
methods  or  discretized  numerical  methods.  In  principle,  continuum 
theories  provide  analytical  solutions  which  are  valid  for  any  choice  of 
layer  thicknesses.  However,  numerical  difficulties  are  encountered  in 
the  numerical  evaluation  of  complete  solutions  especially  for  damped 
systems.  The  discretized  methods  are  based  on  finite  element  techniques 
and  offer  easier  and  more  convenient  numerical  solutions  for  viscoelastic 
layered  systems.  However  the  accuracy  of  the  solution  is  affected  by  the 
choice  of  discretization  scheme,  especially  in  the  high  frequency  range. 

In  this  chapter,  the  general  theory  and  characteristics  of  Rayleigh 
waves  in  homogeneous  elastic  and  viscoelastic  half  spaces  are  briefly 
discussed  in  Section  4.2.  In  Section  4.3,  the  discretized  method  of 
treating  Rayleigh  waves  in  a  layered  system  with  a  rigid  base  is  briefly 
reviewed,  and  several  methods  are  proposed  to  extend  this  method  for  the 
approximate  solution  of  layered  systems  resting  on  a  viscoelastic  half 
space.  The  selection  of  the  fundamental  Rayleigh  mode  is  also  explained 
in  this  section.  The  discretized  method  for  Love  waves  in  a  layered 
system  with  a  rigid  base  is  briefly  presented  in  Section  4.4. 

Several  published  solutions  were  used  for  verification  of  results 
from  the  computer  programs  SITE  and  LOVE  which  were  developed  as  part 
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of  the  research  described  herein  for  site  response  analysis  by  the 
discretized  method.  The  comparative  study,  presented  in  Section  4.5, 
confirms  that  the  discretized  method  can  be  used  for  analyses  of 
surface  wave  propagation  in  a  viscoelastic  layered  half  space. 

Additional  analyses  of  the  seismic  response  of  multi-layered  sites 
excited  by  Rayleigh  waves  will  be  presented  in  Chapter  5. 

4.2  Rayleigh  Wave  in  a  Viscoelastic  Half  Space 

It  is  assumed  that  a  simple  harmonic  wave  train  with  motions  in 
the  xz-plane  only  travels  in  the  x-direction  such  that  the  motion  is 
independent  of  the  y-coordinate  and  that  the  amplitude  of  this  motion 
decreases  asymptotically  with  the  distance  z  form  the  free  surface. 
Haves  satisfying  these  conditions  are  called  Rayleigh  waves  and  were 
first  studied  by  Rayleigh  (1885).  A  solution  corresponding  to  this 
definition  may  be  derived  from  the  general  equations  of  motion  for 
two-dimensional  waves  presented  in  Chapter  3.  Since  the  boundary 
condition  at  z»<*>  requires  that  the  wave  potentials  approach  zero  as  z 
approaches  infinity,  the  solutions  for  Rayleigh  waves  can  be  written: 

<b  *  A  e*qz  e"ikx  (4.1a) 

Y  »  B  e"sz  e"ikx  (4.1b) 

where 

2  2  1/2 

q  *  (k  -  k  p)  ;  taking  the  principal  value, 

2  2  1/2 

s  *  (k  -  k  )  ;  taking  the  principal  value, 

s 

and  A  and  B  are  unknown  complex  constants. 

Equation  (4.1),  in  connection  with  Eq.  (3.2)  defines  the  form  of 
the  displacement  and  stress  field  in  the  half  space.  By  introducing 
the  boundary  condition  that  all  stresses  must  vanish  on  the  plane,  z*0, 
the  following  relationship  between  the  wave  numbers  involved  can  be 


obtained: 


By  squaring  each  side  of  Bq.  (4.2)  and  using  the  relationships 

2  *  *  2 
X  -  <kaA)  -  <W 

*  -  <vy2  -  <yv2 

the  following  equation  is  obtained: 

X3  -  8  X2  ♦  (24  -  16  Y)  X  +  16  (Y  -  1)  -  0  (4.3) 

This  is  of  exactly  the  sane  form  as  the  equation  developed  by  Rayleigh 
for  an  undamped  elastic  medium  except  that  real  velocities  are  replaced 
by  complex  velocities  for  a  damped  system.  This  equation  is  a  cubic 
polynomial  with  coefficients  in  the  complex  plane.  It  has  three 
complex  roots,  X,  which  may  not  be  distinct.  The  root  which  satisfies 
the  original  unsquared  equation,  Bq.  (4.2),  provides  the  fundamental 
solution  for  the  surface  wave,  Borcherdt  (1971). 

The  solution  may  be  found  by  Cardan  techniques,  as  shown  by  Hall 

(1964).  The  solution  was  also  carried  out  by  Borcherdt  (1971),  who 
.2  .2 

showed  that  if  V  /V  is  a  root  of  the  complex  Rayleigh  equation  such 

*2.2  *2  .2 

that  0  < j V  j  /jv  |  <  1  ,  then  V  /V  also  satisfies  Bq.  (4.2). 

r  s  r  s 

This  restriction  on  the  roots  of  Eq.  (4.3)  is  the  same  for  the  damped 
and  the  undamped  case  and  is  used  for  the  selection  of  the  root 
corresponding  to  the  Rayleigh  wave.  In  the  undamped  case,  the  solution 
will  consist  of  three  real  roots  for  Poisson’s  ratio  less  than  0.26, 
and  there  will  be  one  real  root  plus  a  conjugate  pair  of  complex  roots 
for  Poisson's  ratio  larger  than  0.26. 
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Once  the  ratio,  X,  between  the  Rayleigh  and  shear  wave  velocities 
has  been  determined,  the  normalized  displacements,  mode  shape,  of  the 
Rayleigh  wave  may  be  obtained  by  substituting  Eqs.  (4.1)  into  Eqs. 

(3.2). 

The  mode  shape  or  normalized  amplitude  distributions  as  a  function 
of  dimensionless  depth  for  four  different  values  of  Poisson's  ratio  are 
shown  in  Fig.  4.1  for  the  case  of  undamped  elastic  media.  The 
dimensionless  depth  is  defined  as  the  actual  depth  divided  by  the  wave 
length  which  is  inversely  proportional  to  frequencies.  The  particle 
motion  at  different  depths  is  shown  in  the  same  plot.  Since  the 
horizontal  and  vertical  components  of  Rayleigh  waves  are  out  of  phase 
by  the  angle  v/2,  the  trajectories  of  the  particle  motions  are 
ellipses.  The  magnitude  and  direction  of  the  elliptical  motion  is 
dependent  on  depth.  The  following  character itic  may  be  summarized  from 
the  plot  shown. 

1.  The  horizontal  amplitude  decays  rapidly  with  depth  near  the 
surface  and  becomes  zero  at  a  depth  of  approximately  one  fifth  of  the 
wave  length.  The  maximum  negative  amplitude  occurs  at  a  depth  of 
approximately  two  fifths  to  one  half  of  the  wave  length  and  then 
gradually  decays  to  zero. 

2.  The  vertical  amplitude  first  increases  slightly  with  depth  and 
reaches  its  maximum  value  at  a  depth  of  0.05  to  0.15  times  the  wave 
length  below  which  it  decays  rapidly  to  zero,  except  for  materials  with 
zero  Poisson's  ratio  for  which  the  maximum  vertical  displacement  occurs 
at  the  surface. 

3.  The  major  horizontal  and  vertical  disturbances  associated  with 
Rayleigh  waves  are  concentrated  within  one  wave  length  from  the  surface. 
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Fig.  4.1  Normalized  Amplitudes  Versus  Dimensionless  Depth  for  R-Waves 


4.  The  amplitudes  at  depth  decrease  as  Poisson's  ratio  decreases. 

5.  For  wave  propagation  in  the  positive  x-direction,  the 
elliptical  particle  motion  at  the  surface  proceeds  counterclockwise. 

At  the  depth  z  *  0.2  X  where  the  horizontal  amplitude  changes  sign, 
the  direction  of  rotation  reverses.  The  minor  axes  of  the  ellipses  are 
perpendicular  to  the  free  surface  of  the  half  space;  i.e.  the  vertical 
motion  is  stronger  than  the  horizontal  motion  at  all  depths. 

The  distribution  of  the  stress  components  with  depth  is  shown  in 
Fig.  4.2.  The  curves  were  calculated  for  Poisson's  ratio  equal  to  0.25 
(dashed  line)  and  for  Poisson's  ratio  equal  to  0.34  (solid  line).  It 
is  apparent  from  the  plot  that  o  changes  sign  at  z«0.25Xr,  whereas 
<5zz  and  tzz  reach  their  maxima  at  approximately  z/Xf  «  0.3  and  then 
falls  off  exponentially  with  depth. 

Rayleigh  waves  in  a  viscoelastic  half  space  have  been  studied  by 
Borcherdt (1971) .  Computed  mode-shapes  for  the  special  case  of 
Poisson's  ratio  equal  to  0.35  are  shown  in  Fig.  4.3.  It  may  be  seen 
from  this  figure  that  the  effect  of  damping  on  the  distribution  of 
motions  with  depth  is  insignificant  for  the  magnitudes  of  damping 
usually  encountered  in  practice.  Another  small  effect  of  damping  is  a 
slight  tilt  of  the  elliptical  orbits  of  particle  motion.  The  major 
effect  of  damping  is  that  the  waves  will  decay  exponentially  as  they 
propagate  in  the  x-directions.  The  decay  factor  is  approximately 
exp  (-2x8)  per  wavelength  in  the  x-direction. 

4.3  Rayleigh  Waves  in  a  Layered  System 

Rayleigh  wave  propagation  in  layered  media  is  of  great  interest  to 
seismologists  and  has  been  treated  in  standard  textbooks.  A  complete 
theoretical  summary  may  be  found  in  Ewing  et  al.  (1957)  and  a  brief 
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4.2  Rayleigh  Have  Stress  fig.  4.3  Normalized  Amplitude  Versus  Depth  for  Haterlals  with 
Distribution  for  vO  25  Different  Values  of  Shear  Have  Damping  Ratios 0 
and  v-0.34  (after  Viktorov.  <after  Borcherdt,  1971) 


description  of  geophysics  applications  has  been  given  by  Grant  and  West 
(1964).  The  equation  of  motion  for  a  system  of  elastic  layers 
overlying  an  elastic  half  space  was  formulated  in  matrix  algebra  by 
Thomson  (1950)  and  improved  by  Haskell  (1953).  An  equivalent 
formulation  which  makes  calculations  possible  for  higher  frequencies 
was  presented  by  Knopoff  (1964),  Dunkin  (1965)  and  Thrower  (1965).  A 
modification  of  the  matrix  formulation  to  give  faster  machine 
computations  for  modal  solutions  to  a  layered  half  space  was  provided 
by  Watson  (1970) .  Formulations  which  include  damping  were  later 
presented  by  Boncheva  (1977)  and  Silva  (1978).  The  fundamental 
approach  for  all  of  these  methods  is  based  on  continuum  theory  which 
eventually  leads  to  a  complicated  nonlinear  eigenvalue  problem. 

Solution  of  this  problem  involves  serious  numerical  difficulties  in 
many  cases  and  is  complicated  by  the  fact  that  in  layered  systems 
infinitely  many  Rayleigh  waves  (modes)  can  exist  simultaneously. 

A  lumped  mass  finite  element  formulation  tor  a  multi-layered 
system  with  rigid  base  was  developed  by  Lysmer  (1969).  This  method 
leads  to  a  simple  quadratic  eigenvalue  problem  which  can  be  transformed 
to  a  linear  eigenvalue  problem  of  the  double  dimension.  This  problem 
can  be  solved  completely  by  standard  techniques.  Lysmer 's  method  was 
extended  by  Waas  (1972)  to  include  a  consistent  mass  formulation  and 
Love  waves.  This  method,  which  will  be  briefly  reviewed  below,  is  the 
basic  numerical  method  employed  in  the  research  described  herein.  As 
part  of  this  research  several  methods  will  be  introduced  which 
facilitate  the  use  of  the  method  for  cases  involving  a  layered  system 
supported  on  a  viscoelastic  half  space.  These  methods  will  be 
discussed  in  Section  4.3.2.  Methods  for  identifying  the  fundamental 
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node  among  the  many  Rayleigh  waves  which  can  exist  in  a  layered  system 
will  be  presented  in  Section  4.3.3. 

4.3.1  Layered  System  With  Rigid  Base 

Consider  the  semi-infinite  layered  system  shown  in  Fig.  4.4a.  All 
motions  occur  in  the  xz-plane  and  any  point  in  the  system  has  two 
degrees  of  freedom.  The  layered  system  is  treated  as  a  continuum  in 
the  horizontal  direction  but  is  discretized  in  the  vertical  direction 
by  assuming  that  the  displacement  is  continuous  at  each  interface  and 
varies  linearly  within  each  layer.  As  shown  by  Wass  (1972),  the 
equation  of  motion  for  an  n-layer  system  may  be  written  as: 

([A]  k2  +  i(p]k  +  CG]  -  (d2Cm3)  M  -  {0 }  (4.6) 

In  this  equation,  jvj  is  a  vector  containing  the  2n  layer 
interface  displacements  and  CA] ,  M,  CG],  and  [M]  are  the  2n  by  2n 
matrices,  assembled  by  addition  of  layer  submatrices  as  shown  in  Fig. 
3.21  for  the  case  of  body  waves  except  that  the  last  two  rows  and 
columns  of  each  total  matrix  are  not  used  because  of  the  assumption  of 
a  rigid  base  for  which  the  displacements  are  zero.  The  submatrices  [A] 
and  CG]  for  each  layer  may  be  expressed  in  terms  of  complex  shear 
moduli,  constrained  moduli,  and  the  layer  thicknesses  as  shorn  in  Bqs. 
(3.51)  and  (3.53)  and  the  submatrix  [B]  is  redefined  as  follows: 
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Fig.  4.4b  Computation  Model  for  Layered  Soil  System  over  Half -Space 


Bach  sufamatrix  of  [m]  may  be  expressed  in  terms  of  consistent  mass, 
lumped  mass  or  a  combination  of  both  as  shown  in  Eq.  (3.56).  For  a 
given  frequency  u),  it  is  convenient  to  introduce  a  matrix  [c]  such  that 
[C]  •  [gJ  -  w2  [m].  This  reduces  Eq.  (4.6)  to: 

([a]  k2  +  i|Y]k  +  [c])  {v}  -  {0}  (4.8) 

This  is  a  quadratic  eigenvalue  problem  which  has  a  solution  {v}  if,  and 
only  if,  the  determinant  of  the  coefficient  matrix  vanishes.  Bence, 
for  any  given  frequency  the  secular  equation: 

|[A]k2  +  i  [b]  k  +  [c]!  -  0  (4.9) 

defines  the  possible  wave  numbers  for  Rayleigh  waves  in  the  layered 
system.  A  numerical  technique  for  finding  the  eigenvalues  and  the 
corresponding  eigenvectors  in  Eq.  (4.8)  has  been  presented  by  Waas 
(1972).  It  can  be  shown  that  this  equation  gives  4n  eigenvalues, 
kg,  s  •  1, 2, . . . 4n  and  their  corresponding  eigenvectors,  {v}s, 
s  *  1,2, .. . . 4n. 

The  case,  k  ■  0,  can  occur  only  in  an  undamped  system  and  when  the 
frequency  is  equal  to  one  of  the  natural  frequencies  of  the  layered 
systems.  In  this  case  the  motion  consists  of  vertically  propagating 
P-  or  S-waves.  The  particle  motion  is  vertical  or  horizontal. 

The  case,  k  real,  can  occur  only  in  an  undamped  system.  The 
motion  is  similar  to  Rayleigh's  original  surface  wave  in  that  it 
propagates  with  constant  amplitude  in  the  x-direction  with  the  phase 
velocity  c  •  uj/k.  The  particle  motion  is  generally  elliptical  but 
may  be  linear  at  certain  depths. 

The  case,  k  purely  imaginary,  may  occur  in  both  damped  and 
undamped  systems.  It  corresponds  to  a  motion  which  does  not  propagate 
but  simply  decays  in  the  x-direction.  The  particle  motion  is  linear 
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for  undamped  systems  and  generally  elliptical  for  damped  systems. 
Motions  are  in  phase  at  all  points. 

The  case,  k  ■  k  ♦  ik^ ,  k^  «  0,  may  occur  in  both  daaped  and 
undamped  systems  and  is  the  case  of  most  interest  for  applications. 

The  particle  motions  are  generally  elliptical  and  the  motions  propagate 
in  the  x-direction  with  the  apparent  phase  velocity  •  w/k^  and 
it  decays  as  exp  (k^  x) .  Generally,  kf  and  k^  will  be  of 
opposite  sign,  i.e.  the  motion  decays  in  the  direction  in  which  it 
propagates.  However,  the  unusual  case  may  occur  that  k^  and  k^ 
have  the  same  sign,  i.e.  the  motion  decays  in  the  opposite  direction  of 
the  phase  velocity.  This  is  not  a  contradiction  since  it  can  be  shown 
that  for  these  cases  the  group  velocity,  see  below,  is  negative,  i.e. 
energy  transmission  occurs  in  the  opposite  direction  of  the  phase 
velocity. 

It  can  be  shown,  Waas  (1972),  that  if  k  is  a  solution  to  the 
eigenvalue  problem  in  Eq.  (4.8)  and  { v}  is  the  corresponding 
eigenvector  (mode  shape)  then  -k  is  also  a  solution  and  the 
corresponding  eigenvector  [vj  is  the  adjoint  of  {v};  i.e.  the  vector 
obtained  by  simply  changing  the  sign  of  all  horizontal  components  of 
{v}.  The  physical  significance  of  this  is  that  the  same  motion  can 
propagate  in  both  the  positive  and  negative  x-direction.  In  the 
applications  only  the  modes  which  decay  (propagate  energy)  in  the 
positive  x-direction  are  of  interest.  Thus,  in  a  damped  system,  only 
the  2n  modes  with  k^  <  0  are  of  interest  and  the  general  solution  to 
the  equation  of  motion  may  be  expressed  in  the  form: 
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in  which  is  the  mode  participation  factor  for  mode  s.  Undamped 
systems  are  most  conveniently  handled  by  introducing  a  very  small 
damping  ratio  and  selecting  the  appropriate  modes  by  the  condition,  k^o. 
Stresses  and  Strains 

Once  the  displacements  at  the  interface  of  each  layer  are  determined 
the  strains  at  the  midpoint  of  the  j-th  layer -may  be  easily  obtained.  The 
following  expressions  are  for  the  case  of  a  single  mode;  if  several  modes 
are  considered,  then  the  total  strain  and  stress  can  be  found  by 
superposition.  The  strains  are: 


Cx  “  -ik{<v2j-l  +  V2j+1)/2}  e 

M(v2j+2  -  V2j)/hj)  *’ikX  (4'lla) 

Yxz  *  { (V2 j+1  *  V2j-l,/hj  "  ik(V2j  +  V2j+2>/2l  e'llUt 
Stresses  may  be  obtained  by  the  substitution  of  these  strains  into 
Hooke's  law: 
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Cx 

#  *  * 

*  ■ 

Mj  '  2Gj  Mj  0 

< 

Ez 

* 

Txz 

0  0  G 

Yxz 

w  - 

. 

Group  Velocity 

For  each  mode  Eg.  (4.6)  determines  k  as  a  complex-valued  function 
of  w.  A  plot  of  this  function  is  called  a  spectral  (or  dispersion) 
curve,  see  Figs.  4.5  and  4.6.  The  slope  of  this  curve,  U,  is  called 
the  group  velocity. 


du 

0  «  Tt. 

dk 


■  c  ♦ 


dc 
k  — 
*  dk 


(4.12) 
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where  c  is  the  complex  phase  velocity.  The  group  velocity  controls  the 
amount  and  direction  of  energy  propagation,  see  Eq.  (4.16),  and  can  be 
computed  from  Eq.  (4.6)  by  first  differentiating  by  parts  and  then 
premultiplying  by  {v}T 

{v}T([A]  k2  +  i[B]  k  +  [G]  -  to2  [M])  d{v}  + 

{v}T((2k  [A]  +  i[B])  dk  -  2  wdLJ  [m])  {v}  -  0  (4.13) 

The  first  term  must  vanish  since  upon  transposition  and  remembering 
that  [B]  «  -CB]  and  that  the  other  matrices  are  symmetric,  we  obtain 

{dv}T([X]k2  -  i[B]k  ♦  [G]  -  to2  LM]){v} 
which  is  zero  since  (-k,  {v})  satisfies  Eq.  (4.6).  The  remaining  term 


gives 


du 

n  *  — - 
dk 


{ v}T (2k2 


2kui  {v}T[M]  {v} 


(4.14) 


This  expression  can  be  further  reduced  if  the  mode  shapes  are  normalized 
according  to  Waas  (1972),  i.e. 

{v}T(k2[A]  -  [CD)  {v;  -  2k2 
which  by  Eq.  (4.8),  multiplied  by  {v;  implies  that 

{v}T(2k2[A]  ♦  ik[B] )  {v}  •  2k2  (4.15) 


Thus  the  group  velocity  may  be  computed  by 


iv}V)  {vj 


(4.16) 


which  was  developed  by  Lysmer  and  Drake  (1972)  for  real  modes  but  is 
actually  valid  for  the  general  case. 

Energy  Transmission 

Each  of  the  Rayleigh  modes  in  the  general  motion  defined  by 
Bq.  (4.10)  can  propagate  independently:  this  would  simply  correspond 
to  the  case  when  only  one  of  the  mode  participation  factors,  o  ,  is 


non-zero.  However,  if  more  than  one,  mode  propagating  energy  is 
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transmitted  between  the  nodes.  Thus,  the  amount  of  energy  transmitted 

through  a  vertical  plane,  say  of  x*>0,  depends  on  the  composition  of  the 

mode  participation  factors,  a  ,  s*l,...,  2n. 

s 

This  problem  has  been  solved  by  Waas  (1972)  who  showed  that  the 
rate  of  energy  transmission  per  unit  width  in  the  y-direction  of  x»0  is: 

E  -  |  im  flu}*  [R]  {u}]  (*-17) 

where  u)  is  the  frequency,  Im  indicates  the  imaginary  part  and  {ll }  is 
the  complex  amplitudes  at  x*0.  Hence,  by  Eq.  (4.10) 

2n 

{a}  *  ^as  W,  "  ^V]  M  (4-18) 

s»l 

where  CV1  is  «  matrix  which  contains  the  mode  shapes  in  its  columns  and 

{a}  is  a  vector  containing  the  mode  participation  factors  a  , 

s 

s*lr • • • » 2n» 

The  matrix  [R]  is  the  transmitting  boundary  matrix  developed  by 
Waas  (1972),  i.e. 

[R]  -  i  CA]  [V]  [X]  (Vf1  +  Cd3  (4.19) 

where  [X]  is  a  diagonal  matrix  which  contains  the  waves  numbers,  k  , 

8 

s*l,...,  2n,  on  the  diagonal  and  CD]  is  a  banded  matrix  assembled  from 
the  layer  submatrices 

0  (M*  -  2  G*)  0  “ (M*  -  2  G*) 

*  * 

Gj  0  -Gs  0 

0  (M*  -  2  G*)  0  -<M*  -  2  G*) 

*  * 

Gj  0  -Gj  0 

T 

j  ■  1,...,  n.  this  matrix  is  related  to  C&J  through  [B]  »  [DJ  -  [Dj . 

For  the  special  case  of  an  undamped  layered  system,  it  can  be 
shown  that  real  modes  transmit  the  following  energy  per  time  unit  per 


unit  width  in  the  y-direction 

Es  "  -  2  U  ks  ias!2  (4,21) 

where  the  sign  follows  the  sign  of  the  group  velocity,  see  Eq.  (4.16). 

In  the  undamped  case  complex  and  purely  imaginary  modes  transmit  no 
energy. 

Dispersion  curves  for  the  first  few  real  inodes  of  Rayleigh  waves 
in  an  undamped  homogeneous  layer  over  a  rigid  base  are  shown  in 
dimensionless  form  in  Fig.  4.5.  Since  the  group  velocity  is 
proportional  to  the  slope  of  these  curves  only  points  with  positive 
slope  corresponds  to  waves  which  propagate  energy  in  the  positive 
x-direction.  A  more  complete  picture  of  the  variation  of  the  wave 
numbers  with  frequency  can  be  obtained  by  plotting  both  the  real  and 
imaginary  parts  of  these.  This  has  been  done  in  Fig.  4.6  which 
corresponds  to  the  same  case.  In  this  graph  only  the  spectral  lines 
show  as  full  lines  correspond  to  waves  which  propagate  energy  in  the 
positive  x-direction.  The  broken  lines  correspond  to  motions  which 
decay  in  the  positve  x-direction  but  transmit  no  energy. 

4.3.2  Layered  System  over  Half  Space 

The  theory  sunmarized  above  is  applicable  only  to  a  layered  system 
over  a  fixed  rigid  base.  Actual  sites  are  more  similar  to  a  layered 
system  over  a  viscoelastic  half  space.  This  problem  can  in  principle 
be  overcome  by  using  a  very  deep  model  with  many  sublayers.  However, 
such  an  approach  will  lead  to  very  large  matrices  and  thus  expensive 
calculations.  Hence,  some  effective  method  for  better  approximating 
the  half-space  condition  is  currently  needed.  Three  methods  have  been 
investigated  in  this  research: 
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Viscous  Boundary 

The  first  method  investigated  involved  the  use  of  Lysmer- 
Kuhlemeyer  (1969)  viscous  boundary  condition  (dashpots)  at  the  bottom 
of  the  layered  system.  According  to  this  method  the  existence  of  the 
lower  half  space  can  be  simulated  by  adopting  the  following 
relationship  between  the  amplitudes  of  forces  (stresses)  and  the 
horizontal,  U,  and  vertical,  V,  displacements  at  the  interface  with  the 
half  space 


MPvl-C 


xz 

i 

z 


(4.22) 


where 


[«=] 


=  1J0 


(4.23/ 


*  * 

and  p,  and  corresponds  to  the  mass  density,  S-wave  and  P-wave 
velocities  of  the  half  space,  respectively.  The  above  forces  can  be 
added  to  the  equation  of  motion  for  the  n-layer  systems.  This  results 
in  the  equation 

{[A}  k2  -  i[B]  k  +  <[C]  +  [H])j  { v}  =  {0}  (4.24; 

which  is  similar  to  Eq.  (4.8)  except  that  all  matrices  now  have  the 
dimension  (2n  +  2)  x  (2n  +  2).  [H]  is  the  expanded  matrix 


DO  = 


H 


(4.29) 


For  any  given  frequency,  w,  Eq.  (4.24)  has  the  same  form  as  Eq.  (4.8) 
and  the  eigenvalue  problem  can  be  solved  as  for  that  equation. 

The  method  is  hard  to  justify  theoretically  since  the  Lysmer- 


Kuhlemeyer  boundary  condition  was  originally  developed  for  plane 
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vertically  propagating  waves  and  thus  implies  energy  propagation  in  the 
vertical  direction  which  does  not  occur  within  Rayleigh  waves  in 
undamped  systems.  However,  experience  with  the  method  has  shown  a 
pronounced  improvement  in  the  ability  of  the  discretized  method  to 
simulate  a  uniform  half  space.  Some  results  are  shown  in  Fig.  4.5 
which  shows  (in  full  lines)  dispersion  curves  computed  by  Lysmer  (1969) 
for  a  uniform  layer  over  a  rigid  base.  The  dotted  line  through  the 
origin  is  the  exact  solution  for  a  complete  half  space  and  the  circles 
indicate  points  computed  using  the  above  viscous  boundary  condition. 

The  improvement  is  limited  to  the  lower  modes  however.  Bq.  (4.24) 
predicts  higher  modes  which  do  not  exist  in  a  perfect  half  space. 
Travelling  Wave  Boundary 

When  a  plane  wave  travels  in  the  x-direction  through  a  viscoelastic 
half  space  the  boundary  condition  at  the  surface  can  be  shown  to  be 


xz 


[».] 


c  'V 
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V 

v  J 


f  e 


-ikx 


(4.30) 


where 


N 


(M*A  -  ikX*) 

*_ 

G  C 


* 

M  B 

G*  (6  -  ik)J 


(4.31) 


and  T  ,  a  ,  U  and  V  are  the  stress  and  displacement  amplitudes  at  x-0. 
XZ  z 

The  constants  1 


A  ■  kq  (q  -  s)/(sq  -  k  ) 


B  ■  s  (k2  -  q2)/(sq 


k2> 


C  -  q  (k2  -  s2)/(sq  -  k2) 


D  ■  sk  (s  -  q)/(sq  -  k2) 


(4.32) 


2  2 
2  2  2  *  2  2  2  * 

where  q  *  k  -  uj/V  ,  s  =  k  -  ui/V  are  rather  complicated  function 

P  P 

of  k  and  the  above  boundary  condition  cannot  be  used  directly  in 
connection  with  the  discretized  method  discussed  above.  However,  by 
expanding  the  constants  in  Eq.  (4.32)  into  Maclaurin  series  about  k  *  0, 
the  following  approximation  can  be  obtained 


[HJ  *  [».]  *  1  [Hb]  k  *  [»c] 


(4.33) 


where 


N  *  53 


•  *  * 

(V  +  2v  )  G 
S  p 


*  *  * 

(V  +  2V  )  M 
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w- 


*  *  * 

(i  -  1)  V  +  i V  V 
s  s  p 


(i  -  1)  V  +  iv  V  -  2V 
P  s  p  s 


( 4 . 3  5  i 


and  [H  ]  is  the  matrix  given  by  Eq.  (4.23).  As  for  the  case  of  the 
viscous  boundary  discussed  above  the  expanded  forms  of  the  matrices  in 
Eq.  (4.33)  car.  now  be  added  to  the  expanded  forms  of  the  matrices  [A^, 
[b]  and  [C]  for  the  n-layer  system.  This  results  in  the  following 
equation  of  motion 

{([A]  +  [HJ,  k2  -  i([B]  -  [Hb])k+  ([C]  +  [Hj,}  { v)  =  { 0}  (4.36) 

This  eigenvalue  problem  can,  in  principle,  be  solved  and  should 
lead  to  better  solutions  for  low  values  of  k.  However,  the  matrix 
([A]  +  [Ha]*  is  not  symmetric  and  the  matrix  ([Bj  -  £hJ  )  is  not 
skew-symmetric.  Hence,  the  solution  of  Eq.  (4.36)  is  considerably  more 
difficult  than  the  eigenvalue  problems  stated  by  Eqs.  (4.8)  or  (4.24). 
Further  research  needs  to  be  performed  to  investigate  the  performance 
of  this  approach. 
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Variable  Depth  Method 

As  shown  in  Pig.  4.1  Rayleigh  waves  in  a  uniform  half  space 
attenuate  rapidly  with  depth.  The  table  shown  in  the  same  figure 
indicates  that  for  most  solids  only  small  amplitudes  occur  at  a  depth 
of  1.5  x  the  wavelength,  X^,  and  at  a  depth  of  2  x  the  wavelength  the 
amplitudes  are  insignificant  compared  to  the  surface  amplitudes.  Thus 
it  is  to  be  expected,  that  the  fundamental  Rayleigh  mode  computed  from 
a  discretized  model  with  a  rigid  base  at  a  depth  of  B  «  l.SX^  will  be 
similar  to  the  corresponding  wave  travelling  in  a  complete  half  space. 
This  suspicion  is  confirmed  by  Fig.  4.5  which  shows  that  for  RH  >  2ir; 
i.e.  H  >  X^  the  correct  wave  number  is  computed  from  a  rigid  base 
model.  Further  confirming  evidence  will  be  presented  in  Section  4.6. 

In  a  typical  layered  soil  system  over  a  half  space  it  can  be 

assumed  that  the  shear  wave  velocity  of  the  half  space  will  be  larger 

than  the  velocity  of  the  fundamental  Rayleigh  wave.  Thus  X^  in  the 

above  expression  can  be  safely  replaced  it  Xg  for  the  half  space,  and 

the  half  space  can  be  simulated  by  a  uniform  layer  of  the  thickness 

H  »  1.5X  »  1.5  V  /V  where  v  is  the  frequency  in  Hz.  Subdividing 

s  s 

this  layer  into  9  sublayers  the  element  height  becomes  1/6  Vs/v  which 
as  described  in  Chapter  3  is  sufficiently  small  to  ensure  numerical 
accuracy.  Thus,  no  matter  what  the  frequency  the  underlying  half  space 
can  be  represented  by  9  layers  as  shown  in  Fig.  4.4b.  As  will  be 
discussed  in  Section  4.6,  when  the  half  space  extends  all  the  way  to 
the  surface  or  if  the  surface  layers  are  very  soft  compared  to  the 
half  space  even  better  accuracy  in  the  mode  shape  can  be  obtained  by 
subdividing  the  top  layer  into  two  equal  sublayers  as  indicated  by  the 
dotted  line  marked  "optional”  in  Fig.  4.4b. 
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The  variable  depth  method  is  simple  to  implement  and  will  be  used 
in  the  remaining  part  of  this  dissertation  and  in  the  associated 
computer  programs,  SITE  and  LOVE,  which  at  each  frequency  automatically 
adjust  the  depth  of  the  computational  model  as  indicated  in  Pig.  4.4b. 
The  method  ensures  the  computation  of  good  lower  modes  which,  as  will 
be  discussed  below,  are  the  modes  of  primary  interest  for  the  research 
presented  herein. 

4.3.3  Mode  Selection 

Having  adopted  the  variable  depth  method  the  motions  of  a  layered 
system  over  a  uniform  half  space  are  given  by  Eq.  (4.10)  which 
unfortunately  contains  2n  unknown  mode  participation  factors.  These 
can,  in  principle,  be  determined  by  2n  boundary  conditions;  say  by  a 
set  of  forces  acting  on  the  plane  x  *  0.  This  was  the  method  used  by 
Waas  (1972)  to  compute  the  transmitting  boundary  matrix  00  in 
Eq.  (4.19).  Alternatively,  the  mode  participation  factors  could  be 
computed  from  2n  given  displacement  amplitudes  at  x  •  0.  This  would 
correspond  to  solving  Eq.  (4.18)  for  (a}.  However,  in  the  usual  site 
response  problem,  only  one  control  motion  and  thus  only  one 
displacement  amplitude  is  known  at  each  frequency;  say  the  horizontal 
surface  amplitude  at  x  *  0.  It  is  therefore  not  possible  to  determine 
the  general  motion  from  a  single  control  motion.  A  particular  solution 
can  be  obtained,  however,  if  it  is  assumed  that  only  the  fundamental 
mode  produces  the  motion  at  the  control  point  in  which  case  Eq.  (4.10) 
reduces  to 

{6}  .  a£{v}£  e* (wt-kfX)  (4.37, 

where  a£,  (v}£  and  k£  correspond  to  the  fundamental  mode.  This 
is  a  reasonable  approach  since  it  can  be  expected  that,  if  Rayleigh 
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waves  exist  at  the  site,  the  fundamental  node  is  the  major  contributor 
to  the  notion. 

Fundamental  Mode 

Adopting  this  idea,  the  remaining  problem  is  to  select  the 
fundamental  mode  from  among  the  2n  solutions  to  the  eigenvalue  problem, 
Eq.  (4.8).  For  an  undamped  system  the  fundamental  mode  will  be  among 
the  real  modes  which  are  the  only  modes  which  transmit  energy  in  the 
x-direction.  Among  these  modes  the  fundamental  mode  will  be  the  one 
with  the  largest  wavenumber  k  (shortest  wavelength,  lowest  phase 
velocity) .  This  definition  coincides  with  that  used  by  seismologists. 
For  a  homogeneous  half  space  it  corresponds  to  Rayleigh's  original  wave 
and  the  straight  part  of  the  dispersion  curve  in  Fig.  4.5. 

The  total  number  of  real  modes  which  can  exist  in  a  given  system 
depends  on  the  frequency  of  excitations  and  the  natural  frequencies  of 
the  system.  The  latter,  which  correspond  to  vertical  wave  propagation 
between  the  free  surface  and  the  rigid  base,  can  be  determined  by 
solving  the  eigenvalue  problem  in  u  which  results  from  setting  k  •  0  in 
Eq.  (4.9).  For  a  homogeneous  layer  over  a  rigid  base  these  frequencies 
can  be  read  off  at  k  »  0  in  Fig.  4.5.  In  general,  at  any  given 
frequency,  u>,  the  number  of  real  modes  which  propagate  energy  in  the 
positive  x-direction  will  be  equal  to  the  number  of  natural  frequencies 
below  this  frequency.  Thus,  for  the  special  case  of  a  homogeneous 
layer  over  a  rigid  base  with  H  *  1. (uH/V#  »  9.425)  five  real 
modes  will  exist.  This  can  be  seen  by  counting  the  number  of 
dispersion  curves  which  intersect  the  horizontal  line  •  9.475 

in  Fig.  4.5.  This  observation  implies  that  real  modes  will  always 
exist  when  the  variable  depth  method  is  used.  Since  the  same  line 
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intersects  the  straight  part  o£  curve  in  Pig.  4.5,  a  good 

fundamental  mode  will  always  be  computed  by  this  method. 

Actual  sites  cannot  be  properly  modeled  by  undamped  materials. 

Thus  the  cases  of  major  interest  involve  damped  layered  systems.  For 
such  systems  the  definition  and  selection  of  the  fundamental  mode  is  a 
much  more  complicated  matter.  As  discussed  in  Section  4.4  all  modes 
are  complex,  i.e.  they  have  complex  wave  numbers  with  negative 
imaginary  parts  corresponding  to  decay  in  the  direction  of  wave 
propagation.  One  can  therefore  not  make  a  simple  statement  to  the 
effect  that  the  fundamental  mode  is  the  real  mode  with  the  largest  wave 
number.  On  the  other  hand  the  introduction  of  damping  does  have  some 
simplifying  effect  on  the  computations.  For  example,  as  discussed  in 
Section  4.4,  it  does  simplify  the  selection  of  the  2n  appropriate  wave 
numbers  (the  ones  with  negative  imaginary  parts).  Also,  damping  tends 
to  eliminate  singularities  in  transfer  functions  and  dispersion 
curves.  For  example,  all  of  the  nondifferential  points  on  the 
dispersion  curves  shown  in  Fig.  4.6  will  become  differentiable  points. 

With  the  magnitude  of  damping  which  has  to  be  introduced  in 
practical  problems  the  changes  in  numerical  values  of  wave  numbers, 
mode  shapes,  etc.  from  the  corresponding  values  obtained  by  undamped 
analysis  are  not  large.  This  was  already  indicated  by  the  results 
presented  in  Fig.  4.3  which  show  that  mode  shapes  are  virtually 
unchanged  by  the  presence  of  damping.  As  a  further  illustration  of 
this  point  a  450  ft  homogenous  layer  over  a  rigid  base  was  subdivided 
into  15  sublayers  and  analyzed  by  program  SITE  at  different  frequencies 
using  four  different  damping  values  (6  ■  04,  0.0014,  34  and  104). 
Computed  wave  numbers  are  shown  in  Table  4.1  which  also  show  the  values 


for  the  corresponding  undamped  half  space.  At  each  frequency  the 
number  of  modes  included  corresponds  to  the  number  of  real  modes  for 
the  undamped  systems,  as  can  be  seen  from  the  natural  undamped 
frequencies  of  the  system  shown  below  the  table.  The  number  of  real 
modes  is  equal  to  the  number  of  natural  frequencies  below  the  frequency 
of  excitation  as  discussed  above. 

The  wave  numbers  shown  for  the  damped  cases  are  those 
corresponding  to  low  attenuation.  These  were  selected  by  first 
ordering  all  of  the  modes  in  order  of  the  least  magnitude  of  the 
imaginary  part  of  the  wave  number  and  then  selecting  as  many  as 
indicated  by  the  undamped  case.  The  argument  for  this  procedure  is 
that  if  it  is  true  that  damping  has  a  small  effect  on  the  wave  modes 
then  the  damped  mode  corresponding  to  the  fundamental  mode  in  the 
undamped  case  (which  does  not  decay)  should  be  among  the  low 
attenuation  modes  for  the  damped  case.  That  this  is  indeed  so  can  be 
seen  by  comparison  of  the  wave  numbers  marked  by  an  asterisk  in 
Table  4.1. 

It  can  also  be  seen  from  Table  4.1  that  among  the  modes  selected 
the  fundamental  mode  is  the  one  with  the  largest  real  part  of  the  wave 
number.  This  of  course  is  not  surprising  in  view  of  the  above 
definition  of  the  fundamental  mode  for  the  undamped  case. 

In  view  of  the  above  considerations  the  following  procedure,  which 
also  serves  herein,  as  the  definition  of  the  fundamental  mode,  has  been 
adopted  for  selection  of  the  fundamental  mode  in  a  damped  (or  undamped) 
system: 

a)  Compute  the  undamped  natural  frequencies  of  the  system  (by 


setting  k»0  in  Eq.  4.6) 


Table  4.1  Effect  of  Damping  Katlo  on  Wave  Numbers  of  H-waves  In  a  Fixed-base  System 
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b)  Determine  m  ■  the  number  of  natural  modes  be lew  the  frequency 


of  excitation. 

c)  Solve  the  eigenvalue  problem,  Eq.  4.8. 

d)  Sort  the  modes  in  order  of  magnitude  of  the  imaginary  part  of 
the  wave  number. 

e)  Select  from  among  the  first  m  modes  the  one  with  the  largest 
real  part.  This  mode  is  the  fundamental  mode. 

Experience  with  the  above  metho-  *t  a  large  range  of  site 
conditions  has  shown  that  the  above  procedure  and  definition  of  the 
fundamental  mode  leads  to  motions  which  have  all  the  usual 
characteristics  of  fundamental  Rayleigh  waves;  decay  with  depth,  simple 
mode  shape  and  low  phase  and  group  velocity.  For  the  undamped  case  the 
definition  coincides  with  that  used  by  seismologists. 

Least-Decay  Mode 

Several  other  schemes  were  investigated  for  selecting  the 
fundamental  mode  of  Rayleigh  waves.  Among  these  one  method,  herein 
named  the  least-decay  method,  deserves  some  discussion.  Consider  a 
typical  Rayleigh  wave  with  the  wave  number  It  ■  k^  +  iki#  <  0. 
According  to  the  theory  presented  in  Section  4.3  the  wave  has  the 
wave  length  2nAr  and  it  attenuates  as  exp  (kjX)  in  the  direction 
of  wave  propagation.  Hence  the  decay  factor  per  wavelength  is 
exp  (2itki/kr>.  The  ratio  (-ki/kf)  is  thus  a  measure  for  how 
fast  a  wave  decays. 

It  is  to  be  expected  that  the  fundamental  mode,  which  does  not 
decay  in  the  undamped  case,  will  have  a  very  small  attenuation  for 
damped  cases.  The  idea  thus  arose  to  define,  for  damped  cases,  the 
fundamental  mode  as  the  mode  with  the  smallest  value  of  the  ratio 
(-ki/kf).  This  is  the  least-decay  method. 


In  order  to  test  this  method  on  the  case  discussed  above  the  modes 
were  sorted  according  to  increasing  values  of  the  above  ratio.  The 
result  of  this  scheme  is  shown  in  Table  4.2  which  also  shows  values  of 
the  ratio  (“kj/(k  0)).  As  can  be  seen  from  this  table  the  least- 
decay  method  works  perfectly  for  this  case. 

Table  4.2  also  illustrates  the  interesting  fact  that  for  all  modes 

the  ratio  (-k^/fk^))  is  independent  of  the  damping  ratio  and  for 

the  fundamental  mode  this  ratio  is  near  unity.  The  reason  for  this 

follows  from  the  continuum  theory  for  Rayleigh  waves  in  a  homogeneous 

half  space  presented  in  Section  4.4.  Suppose  0  *  0  *  0,  then 

s  p 

simple  substitution  of  complex  wave  velocities  V  *  *  V  (l-i0)  and 

s  s 

V  *  *  V  (l-i0)  into  Eq.  4.2  will  show  that  k  +  ik.  ■  k(l-i0), 

P  P  r  1 

where  k  is  the  wave  number  for  the  undamped  case.  Thus,  by  separation 
of  the  real  and  imaginary  part  and  division,  it  follows  that 
— kf/ (kr6)  =  1.  The  deviation  of  this  ratio  from  unity  in  Table  4.2 
is  therefore  a  measure  of  the  inaccuracy  with  which  the  discretized 
models  represent  the  half  space. 

In  spite  of  the  good  results  which  the  least-decay  method  achieved 
in  this  case  it  was  not  adopted  as  the  method  for  selecting  the 
fundamental  mode.  This  is  so  because  in  cases  involving  a  soft  highly 
damped  layer  over  a  stiff  half  space  with  low  damping  it  does  not 
select  a  mode  which  agrees  with  the  definition  of  the  fundamental  mode 
used  by  seismologists.  Rather,  it  tends  to  select  a  mode  which 
corresponds  to  the  classic  Rayleigh  mode  in  the  half  space  without  the 
surface  layer.  This  mode  may  attenuate  slower  than  the  fundamental 
mode  but  will  have  a  longer  wavelength.  In  fact,  as  will  be  explained 
in  Chapter  5,  this  least-decay  mode  may  be  of  more  interest  to 
engineers  than  the  fundamental  mode. 
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4.4  Love  Waves 


Love  waves  which  are  of  the  form 

Uy  *  U(z)  exp  i(wt-kx)  (4.25) 

with  u  *  u  *  0  can  exist  only  in  layered  systems.  The  simplest 
X  z 

case  involving  a  layer  of  thickness  H  over  an  elastic  half  space  has 
been  studied  by  Love  (1927)  and  Bullen  (1963).  Continuum  methods  for 
the  evaluation  of  Love  wave  mode  shapes  and  wave  numbers  in  undamped 
multi-layered  systems  over  a  half  space  have  been  presented  by  Haskell 
(1953)  and  Ewing  et  al  (1957). 

As  was  the  case  for  Rayleigh  waves  infinitely  many  generalized 
(complex)  Love  modes  can  exist  in  a  layered  system.  However,  at  any 
given  frequency  real  modes,  which  are  the  only  ones  usually  considered 
by  seismologists,  can  exist  only  when  the  material  properties  satisfy 
certain  relations. 

For  multi-layered  viscoelastic  systems  the  mode  shapes  and  wave 
numbers  are  most  conveniently  evaluated  by  discretized  methods  similar 
to  that  used  for  Rayleigh  waves  above.  Such  a  method  has  been  described 
by  Lysmer  and  Waas  (1970)  and  Waas  (1972).  As  was  the  case  for  Rayleigh 
waves  this  method  assumes  linear  variation  of  displacements  within 
layers  and  the  existence  of  a  rigid  base  at  some  finite  depth  which  can 
be  varied  with  frequency  to  ensure  proper  simulation  of  a  half  space. 

For  an  n- layer  system  of  the  type  shown  in  Fig.  4.7  the  method 
leads  to  the  equation  of  motion 

([A]k2  +  [G]  -  a)2  [M] )  { v }  *  {0}  (4.26) 

where  [A],  [G]  and  [M]  are  symmetric  n  x  n  matrices  which  may  be 
assembled  from  sublayer  matrices  as  shown  in  Fig.  3.23  for  the  case  of 
SH-waves,  except  that  the  last  row  and  column  of  each  total  m.  -_rix  are 
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not  used  because  of  the  rigid  base  assumption.  The  submatrices  are 
defined  by  Eqs.  (3.72)  to  (3.75). 

Equation  (4.26)  states  an  eigenvalue  problem  similar  to,  but 
simpler  than,  Eq.  (4.8)  for  Rayleigh  waves.  For  any  given  frequency,  u, 
this  problem  can  be  solved  by  standard  methods  for  the  n  eigenvalues 
k^,  s  *  1,...,  n  and  corresponding  mode  shapes,  { v} s .  Assuming 
a  small  amount  of  damping  each  eigenvalue  will  lead  to  a  pair  of  complex 
wave  numbers  +kg.  since  only  waves  which  decay  in  the  x-direction  are 
of  interest,  the  wave  numbers  with  negative  imaginary  parts  are  retained 
and  the  complete  solution  can  be  written  as  a  linear  combination  of  the 
remaining  modes 


n 

{6}  e1 


(wt-k  x) 
s 


(4.27) 


where  {6}  is  a  vector  containing  the  displacements  at  the  layer 
interfaces  and  a^,  s  *  1,...,  n,  are  unknown  mode  participation 
factors. 

As  was  the  case  for  Rayleigh  waves  these  factors  cannot  be 
determined  from  a  single  control  motion.  However,  by  assuming  that 
only  the  fundamental  mode  (selected  as  described  above  for  Rayleigh 
waves)  exists  Eq.  (4.27)  reduces  to 


{6}  =  0 {v}  ei(u,t-kX) 


(4.28) 


and  the  mode  participation  factor  can  be  determined. 

The  non-vanishing  strain  amplitudes  in  the  jth  layer  are 


l.„  ,  ^ 

2ka  (Vj+l  +  V  6 


(4.29) 


y  «  a 

zy 


■  *  4 
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These  values,  at  x«0,  may  be  used  to  evaluate  the  maximum  shear 
strain  in  each  layer  when  iterating  on  the  soil  properties  according  to 
the  equivalent  linear  method. 

The  group  velocity  of  Love  waves  may  be  computed  from  Eq.  (4.16); 
except  that  {v}T  in  this  formula  should  be  replaced  by  {v}T. 

A  computer  program,  LOVE,  has  been  developed  to  perform  the 
required  computations  not  just  for  a  single  frequency  but  for  transient 
motions  as  explained  in  Chapter  5. 

4.5  Numerical  Examples 

In  order  to  verify  the  Rayleigh  wave  and  Love  wave  solutions 
produced  by  programs  SITE  and  LOVE,  respectively,  these  programs  have 
been  applied  to  several  problems  for  which  exact  solutions  have  been 
published.  Numerical  examples  are  also  presented  which  illustrates  the 
sensitivity  of  the  results  to  variations  in  material  properties  and 
geometrical  conditions. 

4.5.1  Uniform  Half  Space 

The  exact  solution  for  Rayleigh  waves  in  an  undamped  uniform  half 
space  is  readily  available  and  was  presented  in  Section  4.2.  The  exact 
mode  shapes  for  Poisson's  ratio  equal  to  0.25  and  0.45  are  shown  as 
solid  lines  in  Fig.  4.8  which  also  shows  several  fundamental  mode 
solutions  obtained  by  program  SITE. 

The  first  solution,  dashed  line,  for  Poisson's  ratio  equal  to  0.25 
was  obtained  from  a  model  which  consisted  of  nine  sublayers  of  equal 
thickness.  The  depth  to  the  rigid  base  was  1.5  times  the  wavelength  of 
shear  waves.  The  solution  is  generally  good  except  near  the  free 
surface  where  the  piecewise  linear  approximation  is  too  coarse  to  model 
the  curvature  of  the  exact  mode  shape.  This  problem  can  be  overcome  by 
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subdividing  the  top  layer { s)  as  shown  by  the  other  two  solutions  for 
the  case  of  Poisson's  ratio  equal  to  0.25. 

As  Poisson's  ratio  approaches  1/2  the  discretized  method  tends  to 
overestimate  the  vertical  displacements.  As  can  be  seen  from  Fig.  4.8 
the  error  is  about  15%  for  Poisson's  ratio  equal  to  0.45.  This  error 
can  probably  be  reduced  by  further  subdivision  and  the  choice  of  a 
deeper  rigid  base.  However,  the  error  appears  to  be  related  to  the 
well-known  problems  associated  with  the  use  of  the  finite  element 
method  for  plane  strain  problems  with  high  Poisson's  ratio.  In  any 
case  the  error  has  been  judged  acceptable  for  engineering  applications. 
The  phase  velocity  predicted  by  program  SITE  is  in  excellent  agreement 
with  the  exact  value  for  all  the  models  used. 

It  is  clear  from  the  above  results  that  as  far  as  the  fundamental 
mode  is  concerned  a  ten-layer  model  with  the  fixed  base  at  a  depth  of 
1.5  x  the  wavelength  of  shear  waves  provides  an  excellent  approximation 
to  an  elastic  half  space. 

4.5.2  Single  Layer  over  Half  Space 

The  characteristics  of  Rayleigh  waves  in  an  undamped  system 
consisting  of  a  single  layer  over  an  elastic  half  space  have  been  studied 
by  Mooney  and  Bolt  (1966).  The  physical  model  with  the  notation  used 
for  system  properties  is  shown  in  Fig.  4.9.  The  curves  shown  in 
Figs.  4.10-4.12  are  the  solutions  produced  by  Mooney  and  Bolt  for  the 
ratio  between  horizontal  and  vertical  displacements  at  the  ground 
surface,  phase  velocity  and  group  velocity,  respectively.  In  these 
graphs,  T  is  the  period  of  the  motion  and  Bj/®!  (■  Vg/V#')  the  shear  wave 
velocity  contrast  between  the  half  space  and  the  surface  layer.  The 
solutions  correspond  to  the  special  case:  Y  *  Y'  *  0.25,  Y  *  162.5  pcf, 

Y*  -  125  pcf  (i.e.  P/P'  -  1.3). 


formal }zed 
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Discretized  models  were  prepared  for  these  values  of  the  S-wave 

velocity  ratio,  V  /V  •  »  2.0,  3.0  and  4.0.  In  each  case  the  surface 
s  s 

layer  was  modeled  by  18  sublayers  and  the  half  space  was  handled  by  the 
variable  depth  method.  These  models  were  then  analyzed  using  program 
SITE.  The  results  obtained  are  indicated  by  dots  in  Figs.  4.10-4.12. 
They  are  in  excellent  agreement  with  the  exact  solution  and  clearly 
demonstrate  the  adequacy  of  the  discrete  method  and  the  associated 
variable  depth  method  for  layered  systems  over  half  spaces  even  for 
modes  beyond  the  fundamental  mode. 

The  equivalent  problem  in  terms  of  Love  waves  has  been  studied  by 

Stoneley  (1955).  His  model  is  shown  in  Fig.  4.13  which  also  shows  the 

discretization  for  the  surface  layer  used  for  the  corresponding 
discretized  model.  The  half  spaoe  was  modeled  by  the  variable  depth 
method  with  10  sublayers.  In  Fig.  4.14  Stoneley's  exact  dispersion 
curves  for  this  case  are  compared  with  points  obtained  from  the 
discretized  model  using  program  LOVE.  Again,  excellent  agreement  was 
obtained. 

4.5.3  Two  Layers  over  Half  Space 

Stoneley  (1957)  also  studied  the  propagation  of  Rayleigh  waves  in 
an  undamped  system  consisting  of  two  layers  over  a  half  space. 

His  model  for  Rayleigh  waves  is  shown  in  Fig.  4.15  which  also  shows 

the  discretization  used  in  the  corresponding  SITE  model.  The  variable 
depth  method  was  used  to  simulate  the  half  space.  Figure  4.16  shows 
computed  amplitude  ratios  at  the  ground  surface,  phase  velocities  and 
group  velocities  for  the  fundamental  mode.  The  agreement  between  the 
exact  solution  and  the  discretized  solution  produced  by  program  SITE  is 


excellent. 


Crourxl  surface 


Ground  Surface 
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CHAPTER  5 

EXAMPLES  OF  SITE  RESPONSE  ANALYSIS 


5.1  Introduction 

In  this  chapter,  the  steady  state  methods  presented  in  Chapters  3 
and  4  will  be  applied  to  solve  a  number  of  transient  site  response 
problems  of  the  general  type  shown  in  Fig.  5.0.  The  results  presented 
have  been  selected  to  illustrate  the  major  effects  of  horizontal  wave 
propagation  in  the  form  of  fundamental  mode  surface  waves  or  inclined 
body  waves. 

In  all  cases  the  site  is  assumed  to  be  horizontally  layered  and  to 
be  underlain  by  a  homogeneous  half  space.  For  cases  involving  surface 
waves  this  half  space  is  modeled  by  the  variable  depth  method,  see 
Section  4.4.2. 

All  materials  are  assumed  to  be  isotropic  and  viscoelastic. 
Nonlinearitics  are  approximated  by  the  equivalent  linear  method,  see 
Section  2.5.  This  implies  that  the  possibility  of  complete  failure  and 
large  permanent  deformations  of  the  site  are  not  considered. 

The  sites  investigated  include  a  typical  rock  site,  a  cohesionless 
(sand)  site  and  a  typical  alluvial  site  with  high  ground  water  level  and 
alternating  layers  of  sands,  silts  and  clays.  Seme  of  these  sites  are 
related  to  the  soil-structure  interaction  problems  discussed  in 
Chapter  6.  Although  horizontal  wave  propagation  is  the  main  theme  of 
this  chapter,  all  sites  have  been  analyzed,  using  standard  deconvolution 
procedures  (Program  FLUSH) ,  for  the  special  case  of  vertically  incident 
body  waves.  Since  these  are  well-known  procedures,  results  will  be 
presented  without  futher  conrtents  in  the  following  section  whenever  a 
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compar iaon  is  warranted.  Such  results  will  be  identified  in  graphs  by 
the  notation;  "S-wave"  or  "P-wave"  whatever  the  case  may  be. 

5.2  Transient  Motions 

All  of  the  aethods  described  in  the  previous  chapters  have  been 
limited  to  steady  state  harmonic  motions.  The  remaining  part  of  this 
dissertation  will  be  dealing  with  transient  motion  of  finite  duration 
which  better  model  earthquake  motion.  This  transition  is  achieved 
through  the  use  of  Fourier  techniques  which  involve  a  discrete  Fourier 
transform,  complex  transfer  function  and  interpolation  on  the  latter  in 
the  frequency  domain.  This  technique,  known  as  the  complex  response 
method,  has  been  used  extensively  in  recent  years  and  has  previously 
been  described  by  Schnabel  et  al  (1973),  Lysmer  et  al  (1974,  1975)  and 
Idriss  et  al.  (1973). 

5.2.1  The  Fast  Fourier  Transform 

The  basic  input  to  any  seismic  analysis  is  a  digitized  control 
motion,  y(t),  which  will  be  assumed  to  be  given  at  N  (even)  pints  at 
the  uniform  time  interval  At.  Under  these  conditions  the  control 
motion  can  be  written 


V  -  w,t 

y(tl  .  Be  >  I«* 

8*0 


where 


,,  2xs  .  ,  N 

s  N*At  '  •*0'1 .  2 


(5.1) 


(5.2) 


The  differentiable  function  defined  by  Eq.  (5.1)  may  be  thought  of 
as  a  smooth  interpolation  function  between  the  given  points  of  the 


control  motions 
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Bquation  (5.1)  is  a  truncated  Fourier  series  which  implies  that 
the  function  y(t)  is  periodic  with  the  period  T  ■  N*At.  Actual 
earthquakes  are  not  periodic.  However,  this  problem  can  be  handled  by 
adding  a  "quiet  zone"  consisting  of  a  limited  number  of  zeroes  to  the 
given  control  motion,  Schnabel  (1972)}  thus  increasing  N  (and  T) .  If 
the  quiet  zone  is  sufficiently  long  the  strong  motion  occuring  at  the 
beginning  at  each  cycle  will  decay  because  of  material  dampihg  before 
the  beginning  of  the  next  cycle.  Thus  the  response  within  each  cycle 
is  virtually  identical  to  that  of  a  single  earthquake. 

The  complex  coefficients,  Ys,  in  Eq.  (5.1)  can  be  computed  from 
the  given  values,  y^  ■  y  (k*At),  k  *  0,1,...,  N-l,  of  the  control 
motion.  By  choosing  the  length  of  the  quiet  zone  such  that  N  is  a 
power  of  2  this  can  be  done  extremely  efficiently  by  the  Fast  Fourier 
Transform  algorithm  developed  by  Cooley  and  Tukey  (1965).  The  inverse 
version  of  this  algorithm  can  be  used  to  convert  from  frequency  domain 

to  time  domain,  i.e.  to  compute  the  y  values  from  the  Y  values. 

K  3 

In  seismic  applications  it  is  usual  to  neglect  the  first  term  of 
the  sum  in  Eq.  (5.1).  This  is  equivalent  to  assuming  that  the  control 
motion  has  a  zero  mean  value. 

5.2.2  The  Complex  Response  Method 

According  to  the  complex  response  method  the  response  of  any 
linear  system  to  the  real  excitation  defined  by  B?.  (5.1)  can  be 
determined  as  the  real  part  of  the  response  of  the  system  to  the 


oomplex  excitation 


which  siaply  states  that  the  excitation  is  a  finite  sun  of  hamonic 
excitations. 


Using  the  methods  developed  in  the  previous  chapters  the  response 

of  any  point  to  each  of  these  harmonics  can  be  expressed  in  the  form 

U  -  H(u)  )  V  (5.4) 

s  s  s 

where  U  is  a  complex  amplitude  and  H(u  )  are  discrete  values  of  a 
s  s 

smoother  transfer  function. 

Since  siqserposition  is  valid  for  linear  systems  the  real  response 
in  the  time  domain  is 


U(t)  -  Re 


ia}gt 


(5.5) 


s«l 

which  is  similar  to  Bq.(S.l)  and  can  be  evaluated  by  the  inverse  Fast 
Fourier  Transform  algorithm. 

Since  the  number,  N,  of  points  in  the  time  domain  is  typically 
1024  or  2048,  up  to  1024  values  are  needed  for  the  transfer  functions 
H(ojs).  However,  since  these  functions  are  smooth  only  30-40  points 
need  actually  to  be  determined  by  the  rather  complicated  methods 
described  in  the  previous  chapters.  The  intermediate  points  can  be 
obtained  by  a  special  interpolation  technique  in  the  complex  plane, 
Lysmer  et  al  (1975). 

5.3  Linear  Rock  Site 

The  first  site  considered  consists  of  a  50  feet  layer  of 
well-cemented  sandstone  over  harder  bedrock.  Typical  properties  for 
such  a  site  are  shown  in  Table  5.1.  These  properties  were  assumed  to 
be  independent  of  shear  strain  amplitude.  Thus  the  analysis  discussed 
in  this  chapter  is  linear. 


5.3.1  Computational  Model 

In  the  computational  model  the  upper  370  ft  of  the  site  were 
represented  by  13  sublayers  as  indicated  in  Table  5.1.  The  half  space 
below  this  depth  was  handled  by  the  variable  depth  method  as  described 
in  Chapter  4.  This  resulted  in  a  computational  model  with  a  total  of 
23  sublayers.  Details  of  the  discretization  are  shown  in  Fig.  5.3. 

This  discretization  easily  satisfies  the  requirements  discussed  in 
Section  3.4.2  up  to  a  frequency  of  20  Hz  which  was  the  cut-off  frequency 
for  all  analyses  discussed  in  this  chapter. 


Table  5.1 

Properties 

of  Linear 

Rock  Site 

Main 

Thickness 

No 

VS 

vp 

Dampi ng 

Layer 

(ft) 

Sublayers 

(fps) 

(fps) 

Ratio 

1 

40 

5 

3600 

2900 

0.02 

2 

10 

1 

3900 

6100 

0.02 

3 

320 

7 

5600 

6600 

0.02 

half 

space 

varies 

10 

5600 

8700 

0.02 

All  unit 

weights  are 

150  pcf • 

5.3.2  Control  Motion 

The  control  point  is  at  the  ground  surface  at  x»0.  The  horizontal 
control  motion  has  a  maximum  acceleration  equal  to  0.75g.  Its  time 
history  of  acceleration  is  shown  in  Fig.  5.1.  As  can  be  seen  from  the 
upper  part  of  this  figure  the  motion  has  a  broad-band  spectrtan  which 
fits  approximately  an  NRC-type  design  spectrum.  In  the  computations 
discussed  below  a  total  of  4096  points  (At-0.01)  were  used  in  the 
Fast  Fourier  computations,  and  frequencies  above  20  Hz  were  not 
considered.  As  a  result  of  this  low-pass  filtering,  the  computational 
maximum  acceleration  was  0.76  g. 


mrmi 


Fig.  5.1  iha  Control  Motion  and  Its  Response  Spectrua 
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5.3.3  Strain  Compatibility 

The  above  control  motion  is  extremely  strong.  Bence,  the 
appropriativeness  of  a  linear  analysis  for  this  case  might  be 
questioned.  Actually,  the  following  conments  will  confirm  that  for 
practical  purposes  the  linear  approach  is  quite  appropriate. 

While  our  current  knowledge  of  the  behavior  of  rock  at  large 
strain  amplitudes  is  sketchy,  approximate  relationships  between 
effective  dynamic  shear  modulus,  damping  ratio  and  shear  strain 
amplitude  have  been  established.  Typical  relationships  of  this  type 
are  shown  in  Fig.  5.2  which  also  shows  the  effective  shear  strain  range 
computed  for  the  linear  rock  site.  It  may  be  seen  from  this  figure  and 
also  from  the  strain  compatible  properties  shown  in  Fig.  5.3  that  the 
maximum  nonlinear  effects  amount  to  a  15%  reduction  in  the  effective 
shear  modulus  and  a  50%  reduction  in  the  damping  ratio  assumed  for  this 
site.  These  effects  are  within  the  range  of  accuracy  with  which 
engineers  can  currently  determine  these  properties  in  the  field. 

Hence,  it  may  be  argued  that  a  nonlinear  analysis  for  this  case  would 
be  a  purely  academic  exercise  and  that  the  likely  changes  in  the 
results  frcm  those  obtained  by  a  linear  analysis  would  be  small. 
Nevertheless,  an  attempt  to  evaluate  the  maximum  nonlinear  effects  for 
this  site  will  be  made  in  Section  5.4. 

5.3.4  Steady  State  Results 

Before  discussing  the  transient  motion  results  it  is  interesting 
to  study  the  behavior  of  steady  state  fundamental  Rayleigh  waves  on  the 
site. 

The  dispersion  curves  shown  in  Fig.  5.4  are  nearly  constant  in  the 
frequency  range  of  interest.  This  indicates  that  the  site  behaves 
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Fig*  5*3  Analysis  of  Rock  Site 


essentially  like  a  uniform  half  space  with  the  properties  of  the 
bedrock.  This  is  not  surprising  since  the  wavelength  of  the  shortest 
Rayleigh  wave  is  about  five  times  the  thickness  of  the  weathered  crust. 

As  expected  the  mode  shapes  of  the  Rayleigh  waves,  shown  in 
Fig.  5.5,  are  similar  to  those  observed  for  a  half  space.  Their  depth 
of  penetration  is  inversely  proportional  to  frequency  and  motions  below 
a  depth  of  one  wavelength  are  insignificant. 

5.3.5  Transient  Results 

The  computed  variations  of  maximum  accelerations  with  depth  below 
the  control  point  are  shown  in  Fig.  5.3  for  both  the  case  of  pure 
Rayleigh  wave  excitation  and  the  case  of  vertically  incident  shear 
waves.  The  variation  is  typical  for  what  would  be  expected  for  a  half 
space.  The  horizontal  motions  are  lower  than  the  vertical  motions. 
Also,  they  attenuate  faster  with  depth  than  those  determined  from  the 
S-wave  analysis. 

The  variation  of  frequency  content  with  depth  is  illustrated  by 
the  response  spectra  shown  in  Fi".  5.6.  Within  the  upper  part  of  the 
site  the  frequency  distributions  of  the  Rayleigh  wave  motion  do  not 
differ  greatly  from  that  determined  from  S-wave  analysis.  At  greater 
depth  the  vertical  R-wave  components  predominate.  They  are  longer  than 
those  determined  by  S-wave  analysis  especially  in  the  low  frequency 
range. 

As  discussed  in  Chapter  4  the  Rayleigh  wave  field  will  attenuate 
in  the  direction  of  wave  propagation.  For  steady  state  waves  the 
approximate  decay  factor  for  the  linear  rock  site  is  exp  (-2x8) *  0.75 
per  wavelength.  The  effect  over  a  traveling  distance  of  1000  ft  is 
clearly  illustrated  by  the  surface  motion  response  spectra  shown  in 


Depth  Below  Ground  Surface  (ft) 
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Fig.  5*5  Normalized  Mode  Shapes  for  Fuadaaental  Rayleigh  Waves 
-  Linear  Rock  Properties 
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rig.  5.7.  The  rate  of  attenuation  indicated  for  this  site  is  hardly 
important  for  engineering  design.  However,  it  should  be  observed  that 
over  a  distance  of  several  thousands  of  feet  all  the  high  frequency 
components  of  the  original  control  motion  will  vanish.  This  is  a 
strong  indication  that  even  on  rock  sites  high  frequency  surface  waves 
cannot  exist  several  miles  from  the  epicenter  of  an  earthquake. 

5.4  Nonlinear  Rock  Site 

As  discussed  above  linear  analysis  of  competent  rock  sites  is 

probably  appropriate  even  for  very  strong  motions.  However,  in  order 

to  investigate  the  maximum  credible  nonlinear  effects  on  the  above  site 

the  following  modifications  were  made  regarding  the  properties  of  the 

site.  First  it  was  assumed  that  the  upper  15  feet  of  the  site  is 

weathered  to  the  point  where  its  low-strain  seismic  wave  velocities  are 

reduced  to  V  <*1500  fps  and  V  -2900  fps  and,  second,  it  was  assumed 
s  P 

that  the  sandstone,  down  to  a  depth  of  370  ft,  disintegrates  during  the 
presmed  earthquake  to  the  point  where  it  behaves  like  sand,  i.e.  its 
modulus  and  damping  ratio  depends  on  strain  amplitude  as  shorn  by  the 
curves  marked  "Sand”  in  Fig.  5.2.  The  above  considerations  lead  to  the 
computational  model  defined  in  Table  5.2. 


Table  5.2  Low-strain  Properties  of  Nonlinear  Rock  Site 


Main 

Layer 

Thickness 

(ft) 

Sublayers 

(No.) 

(fps) 

<&> 

Damping 

Ratio 

1 

15 

2 

1500* 

2900* 

* 

2 

25 

3 

3600* 

6100* 

* 

3 

10 

1 

3900* 

6600* 

* 

4 

320 

7 

5600* 

8700* 

* 

half 

space 

varies 

10 

5600 

8700 

0.02 

All  unit  weights  are  150  pcf. 


*  Computational  values  vary  according  to  Fig.  5.8. 
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Fiq.  5.7  Attenuation  of  Rayleigh  Wave  Motion  with  Travelling  Distance 
-  Rock  Site 


The  model  easily  satisfies  the  discretization  criteria  described 


in  Chapter  3  up  to  the  20  Hz  cut-off  frequency  of  the  analysis 
performed,  even  when  reduced  strain-compatible  moduli  are  used. 

5.4.1  Strain  Compatibility 

The  above  model  was  analyzed  using  the  equivalent  linear  method 
and  the  same  control  motion  as  for  the  linear  rock  site.  This  lead  to 
the  strain-compatible  properties  shown  in  Fig.  5.8.  It  is  evident  from 
this  figure  that  the  nonlinear  effects  and  thus  the  maximum  stresses 
produced  by  a  Rayleigh  wave  field  in  this  site  are  considerably  larger 
than  those  produced  by  the  corresponding  S-wave  field.  Incidentially, 
the  principal  stress  directions  of  the  two  fields  are  also  completely 
different.  In  an  S-wave  field  maximum  shear  stresses  occur  on 
horizontal  planes  while  in  an  R-wave  field  they  tend  to  occur  on  the 
45°  planes,  at  least  within  depths  of  interest  to  engineers. 

It  might  at  this  point  be  argued  that  the  strain-compatible  soil 
properties  determined  from  the  R-wave  analysis  should  be  used  for  all 
futher  R-wave  calculations.  However,  it  might  also  be  argued  that  from 
a  practical  standpoint  it  makes  more  sense  to  use  the  S-wave  compatible 
properties  for  R-wave  calculations.  This  is  so  because  actual  near¬ 
surface  ground  motions  consist  mainly  of  vertically  or  nearly 
vertically  propagating  body  waves  and  it  is  these  waves,  and  especially 
the  shear  waves,  which  produce  the  major  part  of  the  shear  strains  in 
the  ground.  The  additional  strains  produced  by  a  weak  superimposed 
Rayleigh  wave  field  are  too  small  to  influence  the  choice  of  strain- 
compatible  properties. 

The  latter  approach  has  been  used  to  produce  most  of  the  results 
presented  below.  It  has  the  futher  advantage  that  it  facilitates  the 


superposition  of  different  types  of  wave  fields,  an  operation  which 
would  not  be  valid  if  the  fields  were  determined  from  different  strain- 


compatible  models.  Thus,  in  the  following,  unless  otherwise  mentioned, 
it  may  be  assisted  that  S-wave  strain-compatible  properties  were  used  in 
all  computations  involving  transient  motions. 

5.4.2  Steady  State  Results 

Fundamental  Rayleigh  wave  mode  shapes  at  selected  frequencies  in 
the  range  of  interest  are  shown  in  Fig.  5.9.  The  dependency  of  the  mode 
shapes  on  the  choice  of  rock  properties  indicate  two  major  effects  of 
increasing  nonlinearity: 

•  A  significant  decrease  in  the  ratio  between  the  vertical  and 
horizontal  motions. 

•  A  decrease  in  the  motions  at  depth,  especially  in  the  high 
frequency  range. 

The  classic  half  space  theory  far  Rayleigh  waves  would  predict  the 
second  observation.  It  also  predicts  that  vertical  surface  motions  are 
always  larger  than  the  horizontal  motions:  a  prediction  which  does  not 
agree  with  field  observations.  Thus  the  first  of  the  above 
observations  about  the  effect  of  nonlinearities  may  be  part  of  the 
explanation  for  this  discrepancy.  In  fact,  it  is  only  part  of  the 
explanation.  The  smaller  vertical  motions  observed  in  the  field  can 
also  be  explained  as  an  effect  of  stiffness  contrasts  between  the 
surface  layers  and  the  bedrock.  Nonlinearities  tend  to  increase  this 
contrast  when  strong  motions  occurs.  Bence,  the  two  explanations  are 
closely  interconnected. 

Another  effect  of  nonlinearities  (or  layering  if  one  prefers  that 
explanation)  is  to  increase  the  dispersiveness  of  the  site.  This  can 


be  seen  from  the  dispersion  curves  presented  in  Pig.  5.10.  These  curves 
which  were  computed  using  S-wave  compatible  properties  shew  that  the 
modes  shown  in  Pig.  5.9  propagate  with  completely  different 
velocities.  Generally  the  velocity  decreases  with  frequency  from  a 
high  velocity  corresponding  to  the  velocity  of  the  bedrock  to  a  low 
velocity  corresponding  to  the  velocity  of  the  surface  layer. 

The  rate  of  attenuation  in  the  direction  of  wave  propagation, 
exp  (-k2x) ,  increases  rapidly  with  frequency  and  increasing  magnitude 
of  nonlinearity.  This  can  be  seen  from  the  variation  of  wave  numbers 
shown  in  Pig.  5.11.  The  same  graph  shows  that  dispersion,  which  is 
proportional  to  k^,  increases  with  increasing  nonlinearities. 

5.4.3  Transient  Results 

Transient  results  for  the  variation  of  maximum  accelerations  with 
depth  are  shown  in  Fig.  5.8.  The  horizontal  motions  are  similar  to 
those  computed  for  the  linear  rock  site,  see  Pig.  5.3,  independent  of 
the  choice  of  strain-compatible  rock  properties.  The  vertical  motions 
are  smaller  that  those  computed  for  the  linear  rock  site  and  they  are 
as  expected  smallest  when  R-wave  compatible  properties  are  used.  Thus, 
for  this  site,  it  may  be  considered  conservative  to  use  S-wave 
compatible  properties  for  engineering  computations.  Por  this  reason, 
and  for  the  reasons  given  in  the  previous  section,  S-wave  compatible 
properties  will  be  used  in  all  further  computations. 

The  variation  of  frequency  content  with  depth  is  illustrated  by 
the  response  spectra  shown  in  Pigs.  5.12  to  5.14.  The  spectra  are 
quite  similar  to  those  computed  for  the  linear  rock  site,  Pig.  5.6, 
exoept  that  for  the  nonlinear  rock  site  the  vertical  accelerations 
contain  fewer  high  frequency  components. 


Fundamental  Rayleigh  Hode  Ground  Surface 


Fig.  5.10  Rayleigh  Rave  Dispersion  Curves  for  S-Vave  Compatible  Properties 
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Pig.  5,14  Comparison  of  Horizontal  and  Vortical  Coaponsnt  of  R-Wavt 
Motion  at  Top  of  Haif-spacs  -  Nonlinear  Rock  Sits 


The  attenuation  of  the  Rayleigh  wave  field  in  the  direction  of  wave 
propagation  is  illustrated  by  Figs.  5. IS  and  5.16  which  shows  oomputed 
response  spectra  at  different  distances  from  the  control  point.  The 
results  are  similar  to  those  obtained  for  the  linear  rock  site*  Fig.  5.7, 
except  that  the  attenuation  is  considerably  stronger. 

5.4.4  Conclusions  for  Rock  Sites 

A  rock  site  has  been  analyzed  using  two  extreme  models.  The  first 
model  assumed  linear  behavior  and  the  second  extreme  nonlinear 
behavior.  Except  for  some  minor  differences  the  two  models,  when 
analyzed  for  the  same  control  motion  at  a  surface  control  point,  lead  to 
similar  results  within  depths  and  distances  of  interest  to  foundation 
engineers.  The  major  conclusion  of  the  above  study  must  therefore  be 
that  for  engineering  purposes  nonlinear  effects  need  not  to  be 
considered  in  site  response  analyses  of  rock  sites. 

It  was  also  found  that  within  depths  and  horizontal  distances  from 
the  control  point  of  normal  interest  to  foundation  engineers  the 
horizontal  motions  computed  on  the  basis  of  the  assixaption  of  a 
vertically  incident  S-wave  field  are  very  similar  to  those  computed  from 
a  pure  R-wave  field.  This  similarity  does  not  extend  to  the  phase 
difference  between  distant  points  on  a  horizontal  plane. 

As  opposed  to  vertically  incident  fields,  Rayleigh  wave  fields 
attenuate  in  the  direction  of  wave  porpagation.  However,  for  rock  sites 
this  attenuation  is  too  small  to  be  of  interest  to  engineers.  It  might, 
however,  explain  why  high  frequency  Rayleigh  waves  are  not  observed  by 
seismologists. 

5.5  Cohesionless  Site 

As  an  example  of  a  site  for  which  nonlinear  effects  are  important 
it  was  decided  to  study  a  hypothetical  site  consisting  of  128  ft  of 
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flf .  5.15  Reaponaa  Spactra  of  Horizontal  Component  cf  R-«ava 
Motion  at  Ground  Surface  -  Nonlinear  lock  Site 
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uniformly  dense  dry  sand.  This  site  will  be  analyzed  for  the  effects  of 
Rayleigh  waves,  inclined  body  waves  and  a  combination  of  such. 

The  sand  was  assumed  to  have  the  following  properties: 

Unit  weight  *  125  pcf 
Relative  density  ■  75% 

Poisson's  ratio  “0.3 

Experimental  data  for  the  variation  of  the  modulus  and  damping 

ratio  with  effective  dynamic  shear  strain  amplitude  and  confining 

pressure  for  such  a  material  were  presented  in  Chapter  2,  Pig.  2.3. 

This  data  will  be  used  below  in  connection  with  the  equivalent  linear 

method.  The  bedrock  is  considered  as  a  half  space  with  the  strain- 

independent  properties:  Unit  weight  145  pcf;  Poisson's  ratio,  0.2; 

V  «  4000  fps,  V  *  6532  fps  and  damping  ratio,  2%,  for  both  S-  and 
s  p 

P-waves. 

5.5.1  Harmonic  Rayleigh  Wave 

The  site,  discretized  as  shown  in  Pig.  5.21,  was  first  studied  for 
the  effect  of  a  harmonic  Rayleigh  wave  at  the  frequency  2.5  Hz.  This 
frequency  corresponds  to  the  predominant  frequency  of  the  transient 
control  motion  to  be  discussed  in  Section  5.5.1.  The  wave  was 
normalized  to  produce  a  horizontal  acceleration  amplitude  of  0.25  g  at 
the  ground  surface.  This  corresponds  to  the  acceleration  level  of  the 
transient  control  motion  to  be  used  in  the  next  section. 

Two  solutions  are  presented.  The  first  solution,  shown  in  Pig.  5.17, 
involves  what  many  may  consider  common  engineering  approximations  for 
this  type  of  problem:  The  sand  layer,  which  actually  increases  in 
stiffness  with  depth  was  replaced  by  a  uniform  layer  with  the  oonstant 
shear  wave  velocity,  1148  fps,  and  the  constant  damping  ratio,  5%.  The 
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Fig.  5.17  Harmonic  Rayleigh  Have  and  Vertical  Shear  Have  In  Uniform  Layer  over  Half  Space  -  Frequency  2.5 
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shear  wave  velocity  corresponds  to  the  average  value  within  the  profile 
at  small  strains.  The  damping  ratio,  which  turns  out  not  to  be 
unimportant  for  the  conclusions  to  be  drawn  later,  is  an  engineering 
estimate.  Furthermore,  in  the  analysis  no  attempt  was  made  to  adjust 
the  above  values  for  strain  compatibility,  i.e.  a  linear  analysis  was 
performed. 

From  the  solution  one  may  draw  the  conclusions  that: 

e  Rayleigh  waves  produce  larger  vertical  than  horizontal  motions. 

•  Horizontal  Rayleigh  wave  motions  decrease  faster  with  depth  than 
shear  wave  motions. 

Both  of  these  conclusions  are  in  perfect  agreement  with  classic 
half  space  theories  for  Rayleigh  waves.  The  first  is  rarely  in 
agreement  with  field  observations. 

Now  compare  the  above  conclusions  with  the  results  of  the  more 
complete  solution  shown  in  Fig.  5.18.  In  producing  this  solution  the 
sand  profile  was  modeled  as  a  layered  system,  see  Fig.  5.21,  which 
increased  in  stiffness  according  to  Bq.  (2.1)  and  the  equivalent  linear 
method  was  used  with  S-wave  compatible  properties  to  account  for 
nonlinearities.  This  procedure  actually  underestimates  the  nonlinear 
effects  for  the  R-wave  results,  since,  as  will  be  discussed  in 
Section  5.5.5,  R-waves  cause  larger  shear  stresses  in  the  upper  part  of 
the  profile  than  S-waves. 

It  is  clear  from  the  results  shown  in  Fig.  5.18  that  both  of  the 
above  conclusions  made  on  the  basis  of  the  results  from  the  simplified 
model  are  wrong  for  this  case.  As  they  will  be  for  some  of  the 
transient  cases  to  be  presented  below,  details  of  layering  and  non- 
linearities  must  be  considered  in  site  response  analyses  of  soil 
profiles. 
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Perhaps  the  strongest  effect  is  that  of  nonlinearity  which  in  this 
case  change  the  average  shear  wave  velocity  of  the  sand  layer  from  about 
1148  fps  to  about  772  fps,  a  change  which  will  have  a  pronounced  effect 
on  the  behavior  of  both  R-  and  S-waves. 

5.5.2  Transient  Rayleigh  Waves 

The  above  study  was  repeated  using  a  transient  Rayleigh  wave  field 
defined  by  the  motion  shown  in  Fig.  5.19.  This  motion  was  scaled  to 
0.25  g  and  used  as  the  horizontal  control  motion  at  the  ground  surface. 
1024  points  (At=0.02  sec)  were  used  in  the  Fast  Fourier  Transform 
computations  and  frequencies  above  20  Hz  were  neglected.  Results 
compatible  with  Figs.  5.17  and  5.18  are  shown  in  Figs.  5.20  and  5.21. 
They  generally  confirm  the  conclusion  made  in  the  previous  section. 

The  response  spectra  for  the  motions  at  the  control  point, 

Fig.  5.22,  shows  that,  although  the  maximum  acceleration  of  the 
horizontal  and  vertical  components  are  similar,  the  vertical  component 
contains  higher  frequencies.  Fig.  5.23  shows  the  variation  of  frequency 
with  depth.  As  expected  the  major  effect  appears  to  be  a  reduction  in 
amplitude  approximately  proportional  to  frequency. 

The  average  damping  ratio  in  the  sand  layer  is  about  10%  and  it  is 
therefore  to  be  expected  that  the  attenuation  of  the  Rayleigh  wave  field 
on  the  x-direction  is  very  strong.  This  is  confirmed  by  the  results 
presented  in  Fig.  5.24.  Within  only  250  feet  from  the  control  point  all 
motions  above  5  Hz  have  attenuated  to  insignificant  values.  This 
especially  influences  the  vertical  component  which  already  at  a  distance 
of  about  200  ft  becomes  smaller  than  the  horizontal  component  at  the 
ground  surface.  As  discussed  at  the  end  of  section  5.4  the  high 
attenuation  oomputed  is  strong  evidence  that  fundamental  mode  Rayleigh 
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Fig.  5.20  Site  Response  by  Rayleigh  Waves  and  Vertical  Shear  Waves  -  Transient  Motion 
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Fig.  5-21  Site  Response  by  Rayleigh  Waves  and  Vertical  Shear  Waves 
Transient  Motion  in  Layered  Sand  Site 
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waves  are  not  an  important  contributor  of  near-surface  motions  in 
cohesionless  sites,  at  least  not  at  frequencies  higher  than  about  1  Hz. 
5.5.3  Rayleigh  Wave  Stress  Field 

Another  factor  which  may  impede  the  propagation  of  strong  Rayleigh 
waves  in  cohesionless  sites  is  the  characteristics  of  the  near-surface 
stress  field  produced  by  Rayleigh  waves  as  opposed  to  that  produced  by 
vertically  propagating  shear  wave  fields.  These  characteristics  are: 

•  Shear  wave  fields  produce  no  normal  stresses  on  vertical  and 
horizontal  planes  while  R-wave  fields  produce  very  large 
stresses;  especially  on  the  vertical  plane. 

•  In  a  shear  wave  field  the  maximum  shear  stresses  occur  on 
horizontal  and  vertical  planes  while  for  R-wave  fields  these 
stresses  occur  on  the  45°  planes  near  the  ground  surface  and  on 
the  horizontal  and  vertical  planes  of  depth. 

Both  of  these  characteristics  are  confirmed  by  the  data  presented 
in  Fig.  5.25  which  shows  details  of  the  stress  field  corresponding  to 
the  above  nonlinear  transient  solution  for  cohesionless  site.  Similar 
data  are  shown  in  the  right  hand  parts  of  Figs.  5.17,  5.18,  5.20  and 
5.21.  It  should,  however,  be  observed  that  the  stresses  shown  in 
Figs.  5.17  and  5.20  are  maximum  stresses  normalized  with  respect  to  the 
S-wave  stress  at  the  ground  surface,  while  the  stresses  in  Figs.  5.18 
and  5.21  are  maximisn  normal  stresses  on  the  horizontal  plane.  The  data 
is  confused  by  the  fact  that  S-wave  strain-compatible  properties,  see 
Fig.  5.21,  were  used  in  the  R-wave  analysis.  Hence,  the  stiffness  of 
the  upper  part  of  the  sand  layer  is  not  compatible  with  the  high  maxim  hr 
shear  stresses,  x  (R-wave)  in  Fig.  5.25,  developed  by  the  Rayleigh 


wave  field  near  the  surface.  Nevertheless,  it  is  clear  that  these 
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stresses  cannot  be  sustained  by  the  near  surface  sand  which,  due  to  the 
lew  confining  pressure,  have  very  low  strength.  Thus  a  near  failure 
oondition,  which  cannot  be  handled  well  by  the  equivalent  linear  method, 
will  develop  in  the  top  layers.  This  will  further  impede  the 
propagation  of  strong  Rayleigh  waves.  Even  if  the  top  layer  had  soaie 
strength  due  to  cohesion,  tension  cracks  would  develop  due  to  the  high 
normal  stress  on  the  vertical  plane.  Any  Rayleigh  wave  motion  in  a 
cohesionless  site  must  therefore  be  relatively  weak. 

The  high  normal  stress  on  vertical  planes  produced  by  Rayleigh  wave 
fields  may,  even  for  weak  fields,  induce  high  pressures  on  embedded 
structures.  The  problem  will  be  considered  in  Chapter  6. 

5.5.4  Inclined  Body  Waves 

Since  it  is  unlikely  that  Rayleigh  waves  are  of  importance  for 
cohesionless  sites  the  above  site  was  analyzed,  by  the  method  described 
in  Section  3.5.1,  with  the  same  control  point  and  motion,  for  the 
effects  of  inclined  SV-waves.  As  will  be  shown  such  waves  produce 
essentially  the  same  motions  and  stresses  on  the  site  as  vertically 
propagating  S- waves.  For  this  reason  all  analysis  were  performed 
linearly  using  the  S-wave  compatible  model  shown  in  Fig.  5.21. 

The  fixed  base  complex  natural  frequencies  of  the  18-layer  model 
are  shown  in  Table  5.3.  In  this  table  the  columns  marked  "S-waves"  and 
"F-waves”  corresponds  to  horizontal  and  vertical  modes,  respectively. 

The  inclined  shear  waves  arrive  at  the  base  of  the  sand  profile 
through  the  underlying  viscoelastic  half  space  at  the  incident  angles, 

0,  5,  10  and  20  degrees  from  the  vertical  axis.  Site  transfer  functions 
(defined  as  the  absolute  ratio  between  the  amplitude  in  question  and  the 
horizontal  amplitude  at  the  top  of  bedrock)  for  horizontal  and  vertical 
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Table  5*3  Fixed  Base  Complex  Natural  Frequencies  of  Sand  Site 


1.4803 

0.1597 

2.9387 

4.0958 

0.4354 

8.1312 

6.9444 

0.7051 

13.7863 

9.7657 

0.9834 

19.3872 

12.5171 

1.2536 

24.8494 

15.3489 

1.5282 

30.4712 

18.1037 

1.7904 

35.9402 

20.7286 

2.0332 

41.1513 

23.4591 

2.3297 

46.5719 

25.9610 

2.4938 

51.5388 

28.2083 

2.8073 

56.0003 

30.8749 

3.0326 

61.2941 

32.2445 

3.1404 

64.0131 

34.8634 

3.6592 

69.2123 

35.9375 

3.3149 

71.3446 

37.6328 

3.9796 

74.7102 

39.4701 

3.9002 

78.3577 

41.4184 

2.5956 

82.2255 

0.3171 

0.8644 

1.3999 

1.9522 

2.4886 

3.0339 

3.5543 

4.0365 

4.6251 

4.9509 

5.5732 

6.0204 

6.2344 

7.2644 

6.5809 

3.9004 

7.7427 

5.1530 


Motions  at  different  depths  are  shown  in  Pigs.  5.26  and  5.27, 
respectively.  The  horizontal  transfer  functions  are  for  practical 
purposes  independent  of  angle  of  incidence  within  the  range  investigated. 
So  are  the  shapes  of  the  vertical  transfer  functions.  However,  the 
vertical  components  increase  with  the  angle  of  incidence.  The  horizontal 
transfer  functions  exhibit  peaks  at  the  S-wave  natural  frequencies  of 
the  sand  layer  and  the  vertical  transter  functions  exhibit  peaks  at  the 
P-wave  natural  frequencies,  see  Table  5.3.  Considering  that  the 
approximate  velocity  ratio  between  the  sand  layer  and  bedrock  is  equal 
to  5.5,  these  results  are  in  excellent  agreeaent  with  the  transfer 
functions  presented  in  Chapter  3,  Fig.  3.16,  for  a  similar  profile  with 
the  velocity  ratio  4. 

The  fact  that  the  peaks  of  the  vertical  transfer  functions  tend  to 
occur  at  higher  frequencies  than  the  peaks  for  the  horizontal  transfer 
functions  means  that,  in  nature,  there  will  be  a  tendency  for  vertical 
surface  motions  to  contain  higher  frequencies  than  horizontal  surface 
motions.  This,  incident! ally,  is  in  general  agreement  with  field 
observations. 

A  transient  analysis,  using  the  same  control  motion  as  for  the 
Rayleigh  wave  case  discussed  in  Section  5.2.2  but  with  the  control  point 
at  the  top  at  the  bedrock,  produced  the  response  spectra  shown  in 
Figs.  5.28  and  5.29.  As  might  be  expected  from  the  above  amplification 
study  the  horizontal  motions  are  virtually  independent  of  the  angle  of 
incidence,  see  Fig.  5.28,  while  the  vertical  motions  increase  with 
increasing  angles  of  incidence,  see  Fig.  5.29. 

A  second  transient  analysis  in  which  the  control  point  was  at  the 
ground  surface  gave  similar  results.  Horizontal  response  spectra  from 
this  analysis  are  shown  in  Fig.  5.30. 
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Pig.  5.27  Site  Transfer  Function  on  Vertical  Component  of  SV  Wave 

Notion  -  Sand  Site  Response  to  Inclined  SV  Waves  at  Different 
Angle  of  Incidence 
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Fig.  5.29  Response  Spectra  of  Vertical  Component  of  SV  Wave  Motion 

-  Sand  Site  Response  to  Inclined  SV  Waves  at  Different  Angle 
of  Incidence 
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Fi«.  5*30  Ktaponat  Spectra  of  Hcriaor.tal  Component  for  SV  rfavt  »ct 
la  Sand  site  -  Control  r.otlor.  cn  Ground  Surf  act 
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At  points  below  the  ground  water  level/  whenever  the  P-wave  velocity 
computed  by  the  above  procedure  fell  below  5000  fps  it  was  increased  to 
this  value.  This  modification  was  made  to  account  for  the  fact  that  in 


•  saturated  soil  mass  P-waves  cannot  propagate  slower  than  the  velocity 
of  F-waves  in  water.  Physically  this  means  that  if  the  soil  frame  is 
not  stiff  enough  to  carry  the  P-wave  at  this  velocity  the  wave  will 
travel  through  the  pore  water  rather  than  through  the  soil  frame, 
experimentally,  this  phenomenon  is  often  observed  when  seismic 
refraction  tests  for  P-waves  are  carried  out  on  a  soft  site  with  high 
ground  water  level.  In  fact  xt  is  a  commonly  used  method  to  establish 
the  location  of  the  water  table.  Since  the  ground  water  level  was 
assumed  to  be  located  12  ft  below  the  surface  the  above  procedure  lead 
to  the  strain-compatible  P-wave  velocity  profile  shown  as  a  full  line 
in  the  third  column  of  Pig.  5.34. 

5.4.1  Steady  State  Results 

Steady  state  Rayleigh  wave  computations  were  performed  using  the 
discretized  model  and  the  strain-compatible  soil  properties  shown  in 
Pi9.  5.34.  The  dispersion  curves  for  fundamental-mode  Rayleigh  waves, 
Pig.  5.  37,  show  that  the  site  is  highly  dispersive. 

The  mode  shapes  shown  in  Pig.  5.38  indicates  much  higher  vertical 
than  horizontal  motions.  This  is  probably  true  for  this  site  in  view 
of  the  high  Poisson's  ratio  induced  by  saturation.  However,  as 
discussed  in  connection  with  Pig.  4.8,  the  procedure  used  probably 
overestimates  the  vertical  motions  slightly  for  high  values  of  Poisson's 
ratio.  Another  effect  of  the  high  Poisson's  ratio,  but  not  an  error, 
is  the  tstusual  difference  in  smoothness  between  the  horizontal  and 
vertical  node  shapes  at  higher  frequencies.  Transfer  functions  for  the 


Depth  BeloN  Ground  Surface  (ft) 


Horizontal  Mode  Shape  Vertical  Mode  Shape 


Pig.  508  Normalized  Mode  shape  for  Pundaaental  Rayleigh  Mode 
-  Deep  Alluviax  site 
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horizontal  and  vortical  ocaiponanta  and  differ ant  depths  and  distances 
fraa  tha  aurfaoa  control  point  are  shown  in  Pips.  5.19  and  5.40.  They 
indicate  significant  attenuation  of  all  ooaponente  above  2  Rt  within  a 
few  hundred  feet  fraa  the  control  point. 

5.6.2  Transient  Results 

The  transient  S-wave  analysis  produced  tha  acceleration  profile 
shown  in  Pip.  5.34.  The  steady  Increase  below  a  depth  of  200  feet  is 
probably  due  to  the  fact  that  an  unlikely  control  notion  was  used  for 
this  analysis.  As  shown  by  Rayashi  at  el.  <1971)  and  Reed  et  al. 

(1976)  surface  notions  observed  of  deep  alluvial  sites  do  rot  oontain 
as  a any  high  frequency  components  as  indicated  by  the  control  notion 
spectra  in  Pig.  5.36.  In  fact,  due  to  the  low  velocities  short 
wavelength)  and  high  attenuation  of  shear  waves  in  such  *lf*s  high 
frequency  notions  at  depth  are  highly  attenuated  by  the  tire  they  reach 
the  ground  surface.  As  a  result,  if  a  strong  high- frequent  •  canponent 
is  specified  at  the  ground  surface,  the  deconvolved  notion  at  depth 
beccnes  unrealistically  strong  in  the  high  frequency  range. 

Consequently  the  accelerations  at  depth  becomes  unrealistic  lly  high  ss 
shewn  in  Fig.  5.34. 

The  acceleration  profile  for  the  corresponding  Rayleigh  wave 
analysis  is  shown  in  Pig.  S.41.  The  response  spectre  of  t)  horizontal 
and  vertical  component  of  R-wave  notions  at  three  different  ground 
levels  (at  ground  surface,  at  44  ft  and  332  ft  below  ground  surface) 
are  shown  in  Pigs.  5.42  and  5.43.  Notions  were  oorputed  directly 
under  the  control  point  and  at  a  horizontal  distance  of  500  ft  from 
this  point. 
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UIMIk  (unctions  shown  in  Figs.  5.48  and  5.49.  As  described  before 
It  appears  unlikely  that  surface  waves  of  frequencies  higher  than  2  Bz 
can  actually  exist  on  the  site. 

5.7.2  Transient  Results 

The  transient  analysis  produced  the  variation  with  depth  of  the 
asxiaum  accelerations  and  the  S-wave  compatible  strains  shown  in 
Fig.  5.45.  The  results  are  similar  to  those  found  for  the  deep 
alluvial  site.  The  response  spectra  of  the  horizontal  and  vertical 
component  of  Rayleigh  wave  motions  computed  at  several  different 
traveling  distances  are  shown  in  Figs.  5.50  and  5.51.  Again,  strong 
dacay  characteristics  of  Rayleigh  wave  propagating  in  an  alluvial  site 
Clearly  showa  that  the  high  frequency  Rayleigh  wave  is  not  likely  to 
exist  in  an  alluvial  site. 

5.t  ary  for  All  Sites 

What  has  been  demonstrated  in  this  chapter  is  that  the  numerical 
•ethods  developed  in  Chapters  3  and  4  can  indeed  be  applied  to  actual 
field  problem.  This  demonstration  will  continue  in  the  next  chapter 
where  acwe  of  the  results  obtained  will  be  applied  to  a  number  of  soil- 
structure  interaction  problems. 

Perhaps  even  more  important,  some  of  the  results  obtained  in  this 
chapter  provide  important  information  on  the  likelihood  of  certain 
typee  of  wave  fields  actually  existing  in  nature  and  on  the  relative 
importance  of  horizontal  wave  propagation  to  engineering  projects. 

THIS  is  meat  strikingly  damonstrated  by  the  computed  rapid  attenuation 
at  funPmental  mode  Rayleigh  waves  in  soil  sites  which  seems  to 
preclude  tha  existence  of  high  frequency  Rayleigh  waves  on  such  sites. 
f*«a  observation  does  not  preclude  the  possibility  that  higher-mode 
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Fig*  5*50  Response  Spectra  of  Horizontal  Component  of  B-Wave  Notion 
Travailing  at  Jlfferent  Distance  -  Shallow  Alluvial  Site 


Fi«.  5-51  Response  spectra  of  Vertical  Coaponent  of  R-Wave 
Motion  Travelling  at  Different  Distance  -  shallow 
Alluvial  Site 
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HtfWW a— I  radiation  damping  affect.  This  technique  is  identical 
w  MM  seed  in  the  PLOPH  program,  Lyaaer  (1975),  and  will  not  be 
diSSUSdSd  father.  The  rigid  base  at  a  depth  of  370  feet  has  been  shown 
*9  >  Hearn  U97|),  who  cooperated  in  these  calculations,  to  be 

sudttsgeesly  deep  to  simulate  a  half  space  for  this  model.  It  should 
ee  meet  toned  here  that  the  site  response  analysis  associated  with  this 
enmlytts  msploysd  the  variable  depth  method  for  Rayleigh  waves.  Thus 
•as  as eve  tiftd  base  occurs  only  in  the  interaction  model 
*  »  »  The  Results 

The  results  of  the  CRSAM  analyses  are  presented  in  terms  of  2« 
eveeler etioe  response  spectra  for  the  components  of  motion  at  key  nodes 
»f  file  Structural  model.  In  each  case  the  results  of  the  Rayleigh  wave 
shot  ye id  end  the  combined  body  weve  analysis  are  plotted  together  for 

eeey  rare per  iron. 

e*  ran  te  teen  from  Pig.  (.4  the  motions  computed  for  Point  G  at 
"te  OH#  of  the  bee#  slab  are  nearly  identical  for  the  two  cases 
vn*uh  *t  eat  surprising  since  this  point  is  also  the  control  point  in 
*ae  free  field.  Generally,  the  §  *  P-wave  analysis  is  conservative  for 
lee  t  rental  met  iota  while  the  *-**«  analysis  is  conservative  for 
•<*•**»  *•!  eetiene. 

The  her i rental  sot i one  along  the  length  of  the  base  slab  were  all 
'he  as  et  Point  G.  which  is  not  surprising  since  the  slab  is  very 
net!,  hewer,  rooting  induced  some  variations  in  the  vertical 
nett  one  am  illustrated  in  fig.  t.S.  As  expected,  the  Rayleigh  waves 
ptedhte  slightly  here  rooting  than  the  body  wave  excitation.  The 
*e*ti«et  eetidt  at  the  left  end  of  the  slab  (Point  t)  is  larger  than  at 
the  right  mi  r  Point  n .  This  is  so  because  the  Rayleigh  waves  travel 
fro*  left  re  right  ant  them  attenuate  in  that  direction. 
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S+?  V,'ave  Analysis - 


S+?  Vave  Analysis  - 

R-Wave  Analysis  — — - 


Fig.  6.4  Horizontal  and  Vertical  Response  Spectra  at  Node  C 
(  on  Top  cf  the  Slab  )  -  Rock  Site 
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Fig.  6.5  Spectral  Acceleration  of  Vertical  Motion  alonq  the 
Top  of  the  Foundation  Slab  -  Rock  Site 


1 
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Due  to  the  rocking  of  the  beee  slab,  the  Rayleigh  wave  excitation 
producea  higher  horizontal  Motions  at  higher  point a  in  the  structure. 
This  ia  illustrated  by  the  reaponae  spectra  shown  in  Fig.  6.6.  The 
vertical  aotions  at  the  sane  points  are  shown  in  Fig.  6.7.  These 
notions  appear  to  follow  the  sane  trend  as  the  notion  of  Point  6  on  the 
base  slab,  i.e.  the  Rayleigh  wave  field  gives  higher  reaponae  in  the 
law  frequency  range  (<5  Hz)  and  lower  response  in  the  high  frequency 
range  than  the  body  wave  field.  However,  the  differences  are  not  large. 

In  the  entire  analysis  the  largest  difference  observed  between 
results  obtained  f ran  the  two  seisnic  enviroments  occur ed  at  Point  C 
at  the  top  of  the  internal  atructure.  The  response  spectra  ooeputed 
for  the  notions  at  this  point  are  shown  in  Fig.  6.8.  For  this  one 
point  it  appears  that  at  frequencies  higher  than  8  Bz  the  spectrin  for 
the  horizontal  notion  produced  by  the  Rayleigh  wave  field  is  nearly 
double  as  high  as  the  corresponding  spectrin  for  the  body  wave  field. 

In  view  of  the  above  results  it  appears  that  Rayleigh  wave  notions 
nay  be  critical  far  the  design  of  structures  on  rock.  In  naking  this 
statenent,  it  should,  however,  be  rmaabered  that  no  seisnic 
environment  consists  entirely  of  Rayleigh  waves  and  that  such  waves  may 
not  even  exist  in  the  high  frequency  range  (>4  Hz)  where  the  largest 
differences  were  observed. 

6.4  Structure  on  Sand  Site 

The  analysis  discussed  above  was  repeated  using  the  sane 
structural  nodel  but  this  tine  embedded  into  the  sand  site  studied  in 
Section  S.S.  The  oonputationsl  model  used  is  shown  in  Fig.  6.9. 

6.4.1  Free-Field  Motions 

Contrary  to  the  analysis  of  the  rock  site  previously  studied,  at 
this  site  the  location  of  the  control  point  was  found  to  be  of  crucial 
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?lg.  6.6  Response  Spectra  cf  Horizontal  Motions  at  Node  A  and 
-  Rock  Site 
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Since  the  purpose  of  this  investigation  was  to  determine  the 
behavior  of  the  retaining  walli  only  results  which  relates  to  this 
scope  will  be  presented  below. 

6.5.1  Maximum  Accelerations 

The  maximum  accelerations  along  the  retaining  wall  are  shown  in 
Fig.  6.15  for  the  S-wave  case  and  in  Fig.  6.16  and  6.17  for  the  R-wave 
cases.  In  each  case  the  motions  are  compared  with  the  corresponding 
free-field  motion.  Except  for  the  close-in  Rayleigh  wave  case  the 
interaction  effects  appear  to  be  insignificant  (<10t).  This  was 
generally  true  for  all  points  in  the  structural  system.  In  all  cases 
all  points  of  the  wall  had  essentially  the  same  vertical  acceleration 
which  is  reasonable  for  such  a  stiff  structure. 

The  horizontal  accelerations  induced  by  shear  waves  are  generally 
larger  than  those  induced  by  Rayleigh  waves.  This  observation  does,  as 
will  be  shown  below,  not  mean  that  the  shear  wave  field  is  the  critical 
load  case. 

6.5.2  Shear  Forces  and  Bending  Hoisents 

The  computed  maximum  shear  forces  in  the  retaining  wall  are  shown 
in  Figs.  6. IB  and  6.19.  to  expected,  the  forces  for  the  case  of  the 
close-in  Rayleigh  wave  field  are  somewhat  larger  than  for  the  case  when 
the  control  point  is  located  500  ft  from  the  wall.  Much  more 
significant  is  the  observation  that  the  shear  forces  Induced  by  the 
Rayleigh  wove  field  are  several  times  larger  than  those  induced  by  the 
S-wave  fields.  This  is  so  because,  as  discussed  in  Section  5.5.3,  the 
normal  stresses  an  vertical  planes  are  much  higher  in  Rayleigh  wave 
fields  than  they  are  in  vertically  propagating  fields  (even  vertically 
propagating  P-weves  will  generate  ■sailer  stresses  on  vertical  planes 
than  R-waves  within  normal  depths  of  embedment) . 
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Fig.  6.15  Acceleration  Profile  along  the  Retaining  Vail 
In  Case  of  Shear  Wave  Analysis 
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As  shown  in  Figs.  6.20  snd  6.21  ths  situation  is  siailar  for 
bending  aments  in  ths  wall. 

6.5.3  Design  Considerations 

Ths  high  bending  aments  developed  in  the  retaining  wall  turned 
out  to  be  the  aost  critical  itaa  in  the  design  of  the  raft  system.  The 
practical  solution  to  the  problea  was  to  design  the  wsll  for  a  bending 
ament  which  was  emputed  frm  a  field  consisting  of  1/3  Rayleigh  waves 
and  2/3  vertically  propagating  shear  waves.  This  dscision  was  based  on 
the  argments  presented  in  Section  5.5.3.  according  to  which  it  is 
unlikely  that  a  strong  Rayleigh  wave  field  can  exist  on  the  site.  Bren 
then  the  design  ament  turned  out  to  be  several  tiaes  larger  than  the 
ament  ooaputed  by  say  a  FLUSH  analysis  which  ass  vases  vertically 
propagating  waves. 

More  important  than  the  design  decision  aade,  ths  above  analysis 
illustrates  a  case  for  which  even  a  snail  content  of  Rayleigh  waves  in 
a  control  notion  nay  be  critical,  even  though  such  waves  create  mailer 
accelerations  in  the  structure  then  vertically  propagating  waves. 
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